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ABSTRACT 


The existence of multiple solutions to a general 
Problem of Bolza in the calculus of variations is inves- 
tigated. These multiple stationary solutions are of 
several distinct types, each of which is briefly dis- 
cussed . 

The most important type of multiple stationary 
solutions occurs when a complete set of necessary and 
sufficient conditions have not been applied. A suffi- 
ciency test to eliminate some multiple solutions of 
this type is developed. It is shown that a sufficiency 
test can be broken down into a path sufficiency test and 
an endpoint sufficiency test and that once the first 
necessary conditions of the calculus of variations have 
been applied , the path sufficiency test and the endpoint 
sufficiency test can be applied independently. 

Only the endpoint sufficiency test was investi- 
gated. Analytical application of this condition requires 
the analytical integration of a set of nonlinear differ- 
ential equations subject to mixed boundary conditions. 
Since this is often difficult or impossible, an algo- 
rithm is developed for numerically implementing the end- 
point sufficiency condition. A geodetics problem is 

ix 



solved analytically to illustrate the theory and demon- 
strates that the sufficiency condition is an effective 
computational tool, by eliminating certain classes of non- 
optimal solutions from consideration. 

A survey of numerical methods for solving differ- 
ential equations with mixed boundary conditions arising 
from problems in the calculus of variations is presented. 
The generalized Newton- Raphson method for solving such 
problems is developed in detail and digital computer 
programs for implementing it are included. 

The problem of determining minimum fuel, orbital 
transfers for low-thrust space vehicles is treated as a 
final comprehensive example. Transfers are considered 
which begin in an initial circular orbit and terminate 
in any final elliptic orbit in the same plane. In addi- 
tion the final argument of periapsis, the final true 
anomoly, the final range angle, and the final time are 
all unspecified and considered free. Several examples of 
multiple stationary solutions to a given terminal orbit 
are presented. Digital computer programs for solving 
this problem and for Implementing the endpoint sufficiency 
test are included. The results of a comprehensive study 
of fixed mass, single thrust, orbital transfers from 
Initial circular orbits to a wide range of final elliptic 
orbits are presented. 



CHAPTER 1 


INTRODUCTION 


A renewed interest in the theory of the calculus of 
variations has been witnessed in the past decade due in 
part to the development of sophisticated automatic control 
systems. When the concepts of state and control variables 
from the theory of automatic control systems were incor- 
porated in the theory of the calculus of variations , an 
entirely new practical viewpoint resulted. Once left com- 
pletely to the realm of pure mathematicians, the calculus 
of variations, in its new formulation, speaks a language 
easily understood by engineers. With the potential use of 
the calculus of variations made apparent to a much broader 
spectrum of technology, rapid theoretical development has 
been inevitable. 

Significant advances have been made in casting 
heretofore diverse calculus of variations problems into a 
single generalized formulation, e.g., the Problem of Bolza 
(Bliss, 1946, pp . 187-265). Although not entirely com- 
plete, the formulation of the theory of the calculus of 
variations in control notation used by Vincent and Mason 
(1969) is applicable to most problems of engineering impor- 
tance . 


i 



2 

A limitation to the engineering effectiveness of 
the calculus of variations in its current formulation has 
been the lack of a complete set of sufficiency conditions 
expressed in control notation. Although a fairly complete 
set of sufficiency theorems has been developed for the 
classical problems of the calculus of variations (Bliss , 
1946, pp . 235-265), they are formulated in a manner and 
in a mathematical notation which is difficult to apply to 
most practical problems. For the most general problems 
expressed in control notation, a complete set of suffi- 
ciency conditions is not known. The lack of such a com- 
prehensive engineering formulation has placed heavy 
emphasis on the more easily understood necessary condi- 
tions . 

Trajectories which satisfy all of the necessary 
conditions of the calculus of variations are referred to 
as stationary solutions. Stationary solutions are often 
accepted as optimum trajectories without further investi- 
gation. For most elementary problems this practice is 
adequate. However, stationary solutions are only candi- 
dates for the true optimum. When an investigator obtains 
two or more distinct stationary solutions to his problem, 
this point becomes quite clear. Chapter 2 partially 
resolves the problem of multiple stationary solutions by 
developing a sufficiency condition in modern control 



3 


notation for problems with variable endpoints. This suf- 
ficiency condition is related to classical sufficiency 
conditions . In addition, a numerical algorithm will be 
presented which allows an investigator to apply the end- 
point sufficiency condition in cases where the differen- 
tial equations describing the optimum trajectory cannot 
be integrated analytically. 

However complete a theory may be, it is of little 
use to the engineer unless it leads in a direct fashion 
to solutions of real-world problems. Unfortunately, the 
theory of the calculus of variations inevitably requires 
the integration of a set of nonlinear ordinary differ- 
ential equations with mixed boundary conditions. In 
all but the simplest academic problems, analytic solu- 
tions cannot be obtained. In order for the calculus of 
variations to be a useful engineering tool, a practical 
computational algorithm must be developed. The problem 
of finding numerical solutions to differential equa- 
tions with mixed boundary conditions is the subject of 
Chapter 3. 

Chapter K presents an interesting example problem 
in which bounded control gives rise to multiple station- 
ary solutions in a unique fashion. 

The complex problem of optimum low-thrust orbital 
transfers is investigated in Chapter 5* Several types of 
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multiple stationary solutions are exhibited and the end- 
point sufficiency condition is used to distinguish the 
true optimum from among the candidates. In addition, the 
results of' a comprehensive study of minimum time low- 
thrust orbital transfers is presented. The investigation 
concerns the qualitative aspects of transfers from an 
initial circular orbit to any given final elliptic orbit 
with the angle of transfer and relative argument of 
periapsis unspecified. 

1.1 Historical Background 

The historical development of the calculus of 
variations has been sporadic, and not entirely free of 
controversy. Without going into the details of the 
early development, it is appropriate to mention the out- 
standing contributions and references of historical 
importance . 

In 1696 the Brachistochrone problem was proposed 
by John Bernoulli (Ostwald, 1911, no. 46, p. 3)- He 
sought the path which minimized the time required for a 
mass under the influence of gravity to slide without fric- 
tion from one point to another. In 1744, Euler (1911) 
discovered the characteristic differential equation which 
must be satisfied by optimum trajectories. Lagrange analyt- 
ically formalized the work of Euler and extended the scope 
of problems to two independent variables. While Euler 
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and Lagrange considered necessary conditions for a sta- 
tionary curve ; Legendre (1786); through the use of the 
so-called second variation, developed a further necessary 
condition that distinguished maximal extremals from mini- 
mal extremals. Jacobi (1837) constructed an essential 
modification of the Legendre necessary condition, and 
from it deduced a further test, which set limits on the 
range of the independent variable. Clebsch (1858) and 
Mayer (1868) conducted further investigations connected 
with transformations of the second variation. 

Much of modern theory of the calculus of varia- 
tions is based upon proofs employing strong variations 
which were first developed by Weierstrass. The work of 
Weiers trass, although never formally published, is well 
known largely through the publications of his contempo- 
raries j among them, the works of Bolza (1913) , Kneser 
(1900), Forsyth (1927) and Hilbert (1902) have had the 
most lasting influence. The work of Bolza has been of 
particular influence on current investigators . He formu- 
lated a general problem, now known as the Problem of 
Bolza, which included the problems of Lagrange and Mayer 
as special cases (Bliss, 1946, pp. 189-193)- A consider- 
able interest was taken in the Problem of Bolza, and a 
summary of the works of many investigators is presented 
in the comprehensive book by G. A. Bliss (1946). Recently 



6 

Hestenes (1966) has reformulated the Problem of Bolza from 
the classical dependent variables notation to the modern,, 
control variable , state variable notation. He obtains 
only path sufficiency theorems for problems with fixed 
endpoints and does not discuss the sufficiency condition 
for problems with variable endpoints to be developed in 
Chapter 2 . 

1.2 The Problem of Bolza 

Since much of what is to follow depends upon an 
understanding of the Problem of Bolza as formulated in 
control notation (Vincent and Bruschj 1966, pp. 4 - 5 ) * a 
brief statement of the problem in its simplest form is 
appropriate . 

Among the set of all continuous state functions , 
y ± (t) i = l,2,,.,,nj t Q < t < t f (1.2.1) 

and continuous control variable functions 

u^_(t) k=»l,2,..., men (1.2.2) 

satisfying differential equations and end-conditions of 
the form 

y i = f i( y c j ,u k ,t ) J = ( 1 - 2 . 3 ) 

Mno’nf’VV = 0 

l = 1,2, ...,p < 2n + 2 ( 1 . 2 . 4 ) 

find the set which will minimize a sum of the form: 
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J! = s(y l 0 »y lf .» t o ,t f ) + / L(y 1 (t ) ,u k (t ) ,t)dt . ( 1 . 2 . 5 ) 

J ' t o 


Here it Is assumed that the functions f . , L, g, and 
2 

are of class C . In the above and throughout this presen- 
tation, a dot above a variable will be used to represent 
the derivative of the variable with respect to t, the 
Independent variable. Likewise the subscripts o and f 
will Indicate the evaluation of the variable or expression 
at the initial and final value of t, respectively. For 
the sake of brevity, the range of subscripts i, j, k, and 
l will be as given above and will not be repeated in what 
follows . 

Following the conventional method of Lagrange 
multipliers (Bryson and Ho, 1969)5 minimization of the 
augmented function 


t ~ 

J* = g + + / £ L - A i f i + ^y^at. 

t 

o 


( 1 . 2 . 6 ) 


is considered. In the above equation and throughout the 
presentation, repeated subscripts will be used to imply 
summation. Equation Cl. 2 .6) was obtained by adjoining 
equations (1.2.3) and (1.2. 4) to relation (1.2.5) as fol- 
lows : 

(a) multiplying relations (1.2.3) by the vari- 
ables A^Ct), respectively, integrating from 
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t Q to and by adding the sum of the inte- 
grals to expression (1.2. 5), 

(b) multiplying equations (1.2.4) by the param- 
eters and adding the sum of the products 
to expression (1.25). 

It is convenient to define the following functions: 

<J ( y io’ y if’ t o' t f) = e + UeU (1-2.7) 

H[y 1 (t),x i (t),u :k (t),t] = Xjjfh ~ L (1.2.8) 

The function H is often referred to as the Hamiltonian. 
With these definitions, equation (1.2.6) may be written as 

■fcf 

J* = G + J [-H + X^jdt (1.2.9) 

t o 

By considering small variations in the path and endpoints 
about a nominal path, it can be shown that if the func- 
tions u^.(t) and y^(t) are a solution to the Problem of 
Bolza, then they must satisfy the following necessary con- 
ditions (Hestenes, 1966, pp. 346-351): 

Condition I. There exist continuous multipliers 
X^(t) and Hamiltonian function as defined in equation 
(1.2.8) such that : 

(1) the Euler- Lagrange equations. 


X 


1 




( 1 . 2 . 10 ) 


0 


( 1 . 2 . 11 ) 
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are satisfied at every point along the path and, 
(2) the transversality conditions 


+ H ~ 
eg 0 

0 

(1.2.12) 

6g 

dy. ^io 

J io 

= 0 

(1.2.13) 

6C _ 

eg ' 1 ~ 

0 

(1.2.14) 

i if 

= 0 

(1.2.15) 


are satisfied by the endpoints. 

Condition II. The inequality 
H[y i (t)A 1 (t),uJJ°(t),t] < H[y i (t),X 1 (t),u k (t),t] (1.2.16) 

must be satisfied for all t Q < t < t f and for all non- 

no 

optimal control functions u . This is referred to as the 
necessary condition of Weierstrass. 

Condition III. The k by k matrix 
^2 h s = 1,2, . . . ,k 

^?^t t - 1,2, ... ,k (1.2.17) 

must be negative semi-definite for a minimum. This condi- 
tion is known as the Legendre-Clebsch necessary condition. 

Condition IV. A fourth necessary condition is 
discussed by Bliss (1946, pp. 226-228) in classical depen- 
dent variable notation. He proves that the second order 

-X- 

variation of a sum similar to J must be non-negative 

•X* 

along a. stationary arc, if that arc minimizes J . 



10 


Hestenes (1966, pp. 283-286) verifies this conclusion in 
modern control notation for the Problem of Bolza with 
fixed endpoints. As developed by Hestenes, the fourth 
necessary condition represents a necessary condition on 
the path alone; variations in the endpoints are not con- 
sidered. If Condition III is satisfied. Condition IV is 
usually referred to as the Jacobi Condition. A further 
geometric interpretation of this condition is presented 
in the following section. 

1.3 The Nature of Multiple Stationary Solutions . 

It is helpful in understanding the nature of mul- 
tiple stationary solutions to classify them by the cir- 
cumstances pertinent to their occurrence. The classifi- 
cation is not an apriori one, but rather one based on 
experience . 

1.3-1 Fixed Endpoint Problems and Path Suffi- 
ciency Conditions . Multiple stationary solutions are 
often obtained for problems with fixed endpoints, because 
only the first three necessary conditions have been con- 
sidered. Bliss (1946, p. 235) has shown In dependent 
variable notation that the fourth necessary condition of 
Jacobi, taken together with the first three necessary 
conditions, suitably strengthened, forms a sufficient set 
of conditions for the Problem of Bolza. 
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Although analytically complex, the Jacobi condi- 
tion has a simple geometric interpretation for problems 
with one state variable. In this case the set of all 
solutions forms a one dimensional family of extremal arcs, 
y = y(t,c), which all pass through the initial point as 
shown in Figure (l.l). With each arc is associated a 
particular value of c . If arcs y(t,c) and y(t,c 4- e) 
intersect In the limit as € goes to 0, the point of inter- 
section is called a conjugate point. The locus of such 
intersections is called the discriminant locus, also shown 
in Figure (1.1) . 



Fig. 1.1 A Family of Extremals and the Discriminant 
Locus 
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In terms of this geometry, the Jacobi condition 
requires that an optimal trajectory contain no conjugate 
point. Alternatively, the condition requires that an 
optimal solution may not touch the discriminant locus. 
Figure (1.1) shows that there are two solutions joining 
points 0 and A, one of which touches the discriminant 
locus and is therefore non- optimal. If the fourth neces- 
sary condition of Jacobi is applied in such cases of mul- 
tiple stationary solutions, all but one of the trajector- 
ies will be shown to contain a point conjugate to the 
initial point, thus rendering them non-optimal. 

Examples of this occurrence are many. The 
Brachistochrone problem with fixed endpoints graphically 
illustrates the problem. Consider the problem of a bead 
sliding down a wire under the influence of gravity alone. 
What should the shape of the wire be in order to minimize 
the time of transit between two points in a vertical plane? 
It is well known that the solution curves are cycloids. 
However, as shown in Figure 1.2, there are several differ- 
ent cycloids which satisfy the necessary conditions of 
the calculus of variations . It can be seen that the 
x-axis forms the discriminant locus and that the points 
where solutions 1 and 2 touch the discriminant locus are 
conjugate points. Since solutions 1 and 2 violate the 
Jacobi condition, it Is evident that solution 3 is the 
true opt irnum . 
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Pig. 1.2 Multiple Stationary Solutions for the Brachi- 
stochrone Problem 


In this case it has been possible to distinguish 
the true optimum from the candidates by applying suffi- 
ciency conditions pertaining to the path. Although the 
fourth necessary condition of Jacobi has been exhibited 
in control notation for fixed endpoint problems 
(Hestenes, 1966, pp. 250-286) , it has not been frequently 
used in engineering problems due to complications in 
applying it to cases where the Euler-Lagrange differential 
equations cannot be Integrated analytically. 

1.3-2 Variable Endpoint Problems. Problems with 


variable endpoints require that the endpoints of the tra- 
jectories, as well as the path, be selected in an optimum 
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fashion. Consider the geodesics problem of trying to 
find the minimum distance from the origin to a given 
parabola as shown in Figure (1.3)- It will be shown in 
Chapter 2 that two stationary solutions exist., viz.,, OA 
and OB. Once an endpoint Is determined , the problem 



Fig. 1.3 Multiple Stationary Solutions to a Geo detics 
Problem 


becomes one of fixed endpoints, and the path sufficiency 
conditions previously discussed can then be applied. In 
this case, with the endpoint, A or B, thought of as being 
fixed, it can be shown that both solutions satisfy the 
Jacobi sufficiency conditions regarding path (Bolza, 1961, 
pp. 84-86). Since both paths satisfy sufficiency condi- 
tions, but the endpoints were determined from necessary 



conditions only, it is apparent that a sufficiency condi- 
tion regarding endpoints is needed to distinguish the 
true optimal. Such an endpoint sufficiency condition 
does exist and will he derived in Chapter 2. It will be 
shown that only one of the solutions will represent a 
local minimum with respect to variation of the endpoint 
along the parabola. 

1.3-3 Problems Requiring Path and Endpoint Suffi - 
ciency Conditions . In the last two sections , the neces- 
sity of using path sufficiency conditions and endpoint 
sufficiency conditions was illustrated separately. It is 
not unusual, however, to encounter problems requiring the 
application of both sufficiency conditions to distinguish 
the true optimum from the set of multiple stationary solu- 
tions. To illustrate this situation, reconsider the Brach- 
istrochrone problem where the final endpoint, instead of 
being fixed, is required to be on curve E as shown in Fig- 
ure (1.4). Both trajectories OAB and OB satisfy the end- 
point sufficiency conditions; that is, both solutions 
represent a local minimum with respect to small variations 
of the endpoint along endpoint manifold E. As discussed 
before, point A is a conjugate point, thus violating the 
Jacobi path sufficiency condition. Trajectory OC 
violates the endpoint point sufficiency condition; that 
is, endpoint C represents a local maximum with respect to 
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0 A 



Fig. 1.4 A Problem Requiring both Path and Endpoint 
Sufficiency Conditions 


small variations of the endpoint along E. Thus through 
the use of both sufficiency conditions, OB is selected 
as the true optimum. 

1.3-4 Problems with Periodic Solutions . Consider 
a system of equations (1.2.3) which exhibit periodic oscil- 
lations when no control effort is applied. It is not 
unusual for the optimal controls and adjoint variables of 
such a system to also demonstrate periodic motion with 
the same period. This is especially true if the magnitude 
of the control Is small. By restricting the magnitude of 
the control to sufficiently small values, the deviation 
between the solution from one period to the next can be 
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made as small as desired. Problems with bounded control 
often have controls with such small magnitudes . 

The criterion by which the solution is terminated 
is often periodic for problems exhibiting periodic oscil- 
lations. The terminating or cutoff condition is obtained 
from the transversality conditions (1.2.12) - (1.2.15) by 
eliminating the parameters, to form a single relation- 
ship among the state and adjoint variables. The zeros of 
the cutoff function then represent the terminating condi- 
tion. As an example,, consider the problem of a thrusting 
harmonic oscillator: a mass is connected in parallel by a 
spring and dashpot to an inertial reference. The mass is 
capable of generating a bounded thrust in the upward direc- 
tion. The problem is to get the mass to a specified 
height while minimizing the integral of the thrust with 
respect to time. For simplicity it is assumed that the 
mass is constant. This problem is discussed in detail in 
Chapter 4. For small damping factors and null thrust, 
the state variables, position and velocity, the adjoint 
variables, and the cutoff function all exhibit damped 
periodic oscillations. As shown in Figure (1-5) the cut- 
off condition is satisfied during each period. For suffi- 
ciently small thrust amplitudes the cutoff function will 
deviate only slightly from that generated for null thrust, 
and will be satisfied at several points. 
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Time 


Fig. 1.5 Multiple Solutions Due to a Periodic 
Cutoff Function 


Each time the cutoff condition is satisfied ^ a. 
potential optimal endpoint and a corresponding stationary 
solution is obtained. Thus in periodic systems with weak 
bounded control,, multiple stationary solutions may be 
encountered. Examples of such systems are presented in 
Chapter 4 and Chapter 5- 

1.3-5 Problems with Singular Control . In the cas 
of optimization problems with bounded control variables, 
for example, u . > u > u , one encounters the intrigu- 

Ing problem of the possibility of singular control. In 
cases where the Hamiltonian is linear in a bounded con- 
trol variable,, the control cannot be determined from 
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Euler- Lagrange equation (1.2.11). In this case the 
Hamiltonian can be written as 


H = S(y i> X i .t)u + Q(y i >X i ,t) 


( 1 . 6 . 1 ) 


for scaler control. S Is referred to as the switching 
function. The well-known Maximum Principle for problems 
with bounded control developed by Pontryagin et al . (1962) 
requires that 


u = u when S > 0 

max 

u = u . when S < 0 
min 


( 1 . 6 . 2 ) 


However, for the case when S = 0 over a non-vanishing 
time interval, the Maximum Principle is indeterminant, and 
u may take on intermediate values. This is the case of 
singular control. George Leitmann has pointed out that 
( 1966 , pp. 57-58), "While it is possible in a particular 
problem. . .to rule out the possibility of [singular con- 
trol] , this cannot be done in general." Thus, whenever 
the switching function goes to zero, the control will 
change, i.e., the control will switch to the opposite 
extreme or to singular control. 

To demonstrate the existence of multiple station- 
ary solutions In the case of singular control, examine the 
problem of minimizing 


CO 

J = J x-j 2 dt 


o 
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subject to the constraints: 

X 1 = x 2 + u x l(°) = x iq x- L (oo-)'= 0 

(1.6.4) 

x g = - u x 2 (0) = x 2Q x 2 (co) 0 

|u| <1 

This problem was first discussed by Johnson and Gibson 
(1963). 

The Hamiltonian is 

2 2 ■ 

H = (l 1 ~l 2 )u + l 1 x 2 " X 1 / 2 = ( x 2 + u ) “ ^ 2 " x i / 2 

In this case S = Xj - X 2 * If S is identically zero for 
a’ non- vanishing time interval, singular control exists. 

By taking a suitable number of time derivatives, it can be 
shown that the singular control is given by u = -x-^-Xg, 
and that the singular arcs are two lines x-^(t) =0 and 
x 1 (t) 4- 2x 2 (t) = 0. 

Figure (1.6) shows two possible stationary solu- 
tions to the problem starting at point A, one of which 
has a singular subarc. The first arc AB is the same for 
both solutions. In both solutions u = -1 on arc AB. 
However, at point B, two choices can be made for the 
optimal control. One can continue with u = -1 along arc 
BC and then follow arc CO to the origin with u = + 1. 

This is the so-called "bang-bang” solution. Alternately, 
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at point B one may elect singular control., u = x^, and 
proceed to the origin directly along arc BO. 



Fig. 1.6 Multiple Solutions Arising from Singular 
Control 


Unfortunately, there is no guarantee that the 
solution with singular control is minimizing or that it 
will always enter the optimal solution., even if the possi- 
bility of singular solutions does exist. In this case the 
solution with hang-bang control, arc ABCO, has an index of 
performance almost 12 per cent larger than for the true 
optimal control which uses the singular control arc BO. 

Recently, Kelly, Kopp and Moyer (1967) and Robbins 
(1965) have developed a new necessary condition for test- 
ing the optimality of singular subarcs. Although this 
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necessary condition will certainly eliminate a substan- 
tial number of non-optimal singular solutions., a complete 
set of necessary and sufficient conditions still remain 
to be found. 



CHAPTER 2 

AN ENDPOINT SUFFICIENCY CONDITION 

In Chapter 1 it has been shown that failure to 
apply sufficiency conditions to problems in the calculus 
of variations often results in multiple solutions, some of 
which are non-optimal. The examples have established that 
multiple solutions occur in many types of variational prob- 
lems and for the simplest of problem formulations. In 
addition, a complete set of sufficiency conditions needed 
to select the true optimal has not been formulated in 
modern state-variable, control-variable notation. 

In this chapter an endpoint sufficiency condition 
is developed for variational problems with variable end- 
points. Among the multiple stationary solutions that may 
exist, the condition provides a test for distinguishing 
those stationary solution endpoints which represent a mini- 
mum. By way of proving the endpoint sufficiency condition, 
a novel proof of the transversality necessary condition Is 
also exhibited. The endpoint sufficiency condition is 
related to sufficiency conditions of the classical calculus 
of variations in section 2.4. In section 2.5 the suffi- 
ciency condition is illustrated with an example problem. 
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Finally, a numerical algorithm is developed for applying 
the endpoint sufficiency condition to problems with no 
analytic solution. 


2.1 Functional Relationships for the Problem of Bolza 

Heuristic arguments in Chapter 1 have implied that 
the necessary and sufficient conditions for the Problem of 
Bolza fall into two classes: those pertaining to the path 
and those pertaining to the endpoints. It has been further 
implied that sufficiency conditions pertaining to the end- 
points can be considered independently from those pertain- 
ing to path. Consider the Problem of Bolza as expressed 
in section 1.2, equations (1.2.1) - (1.2.5). In this 
formulation, and for the remainder of this chapter,, the 
controls u^ are assumed to be unbounded functions of 
time. In addition, it is now assumed that the Jacobian 
is not equal to zero 


I dH SH 
3 u. 

8 (u-^, Ug , . . . , u fc ) 


5H 

chT-^, chlg, . . . , 'du k 


£ 0 


( 2 . 1 . 1 ) 


for all points (u^,Ug, . . . ,u k ) In the control space. The 
function H has been previously defined in equation 
( 1 . 2 . 8 ). If equation ( 2 . 1 . 1 ) is valid, the implicit func- 
tion theorem (Buck, 1965, pp. 283-286) assures the 
existence of the k functional relations 


u k = u k[yi(t),x±(t)] 


( 2 . 1 . 2 ) 
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because of the k control variable Euler- Lagrange equa- 
tions (1.2.11). Condition (2.1.1) specifically elimi- 
nates from consideration those systems in which any state 
variable derivative , y as defined in equation (1.2.3) > 
is a linear function of any of the control variables . It 
is therefore assured that the control variable Euler- 
Lagrange equations (1.2.11) will be explicit functions of 

u, . 
k 

A solution to the Problem of Bolza is specified 
by solutions for the state variables yv , as well as the 
control variables u^., as functions of time. A selection 
of the initial time and the final time completes the 
solution. To obtain these solutions, the control variable 
Euler- Lagrange equations are first solved for the control 
variables u^. as functions of the state variables y^ and 
the adjoint variables X^. The control variables in the L 
function in equation (1.2.5) and in the fj_ functions in 
the state variable differential equations (1.2.3) are then 
replaced by the functional relationship for u given in 
equation (2.1.2). The functions f ^ and the L function are 
now explicitly dependent only on the state variables, the 
adjoint variables, and time. 

q = f ± [yj(t); q(yj(t),x ; j(t));t] 

L = L[y-j(t);u k (yj(t) Aj(t));t] 


(2.1.3) 

(2.1.4) 
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With relations (2,1.3) and (2.1.4) substituted into the H 
function, equation (1.2.8), the H function becomes an 
explicit function of the state variables, the adjoint vari- 
ables, and time, alone: 


H 


.o 


H°[y, (t);L (t) ;u, (y . (t),X, ( t ) ) ; t ] 


(2.1.5) 


r'"'!' - '" i v / J "i ' 

Finally, by considering the adjoint variable Euler- Lagrange 
equation (1.2.10) in view of equation (2.1.5)* it is evi- 
dent that a similar functional relationship exists for the 
time derivatives of adjoint variables, 


X, 


p,- [y, (t)*x, (t),t] = - 


8h 


( 2 . 1 . 6 ) 




In summary, once the optimal control is selected, 
the state variable differential equations (1.2.3) and the 
adjoint variable differential equations (1.2.10) comprise 
a set of 2n first order nonlinear differential equations 
in the 2n state and adjoint variables and time. This set 
of differential equations can be integrated in theory, 
yielding 




y i (t*c Y> ) 


1,2, . . . , 2n 


q = q(Uc r ) 


(2.1.7) 

( 2 . 1 . 8 ) 


where the c^'s are constants of integration. The p state 
variable constraints (1.2.4) and the 2n + 2 equations 
representing the transversality necessary conditions 
(1.2.12) - (1.2.15) comprise a set of (2n + p + 2) non- 
linear algebraic equations in the 2n constants c , the p 
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parameters , the initial time t Q , and the final time 
t^. These equations may be solved for the c^'s, t Q and 
t^ . However , since the set of equations is nonlinear, a 
unique solution for these quantities is not guaranteed. 

The initial values (y. . , t ) and the final 

values are specified. Hence, the c^'s may 

be determined as a function of initial and/or final val- 
ues by evaluating equations (2.1.7) and (2.1.8) at either 
the initial or final point. For example, a solution for 
the c 's would be specified by the set (y. ,X. ,t ,t„), 
the set (y if ,X if ,t f ,t Q ) or the set (y lo ,y if , t Q ,t f ) . 

While it is difficult to attach any physical meaning to 
the initial or final values of the Lagrange multipliers, 
the initial and final values of the state variables have 
an immediate physical significance. For this reason, the 
state variable endpoints have been selected to function- 
ally represent the c f constants of integration for the 
rest of this chapter. Thus equations (2.1.7) and (2.1.8) 
can be written as 

h = L^’Uo’yif'Vh) (2.1-9) 

and t Q < t < t^ 

L = (s.i.io) 

By substituting the functional relationships exhibited in 
equations (2.1.9) and (2.1.10) into relations (2.1.2) - 
(2.1.6), it can be seen that the functions u-^,L,f^,H° 


3 



28 


and can all be written as explicit functions of the 

set (t.y. .y.„.t ,t„). These functional relations,, 

\ o f' 

together with that for the function G from equation 
(1.2.7) are summarized for reference below: 

u 


H 


.0 


T o, 


k 

= u k (t,y i 0 ,y lf ,t 0 ,t f ) 

( 2 . 1 . 11 ) 

L 

= L (t,y l 0 ,y lf ,t 0 ,t f ) 

( 2 . 1 . 12 ) 


= qft.yio’yif’k'U) 

(2.1.13) 


= yt-yic.ru-to-tf) 

(2.1.14) 

y io 

,y if 'k’h P x j ^ ,y io ,y if ’h’h ^ 



u k( t<y io’ y if' t o ,t f) ; t] 

(2.1.15) 

G 

= °(>V y io’ y ±f'V t f> 

( 2 . 1 . 16 ) 


So that there will be no confusion as to the meaning of 
the subscripts ,, note that 

(2.1.17) 


y . = y . 


t=t 


o 


X JO ' b 


t=t. 


( 2 . 1 . 18 ) 


o 


y jf s y J 


t=t; 


(2.1.19) 


X jf " X j 


t=t f 


( 2 . 1 . 20 ) 


Using the functional relationships summarized 

•¥: • 

above form the aiigmented function J - J + p t where 

JJ ' JO 
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( 2 . 1 . 21 ) 


-t 


o' f ' 


ft 


+ 


-H°(t,y, .y. f ,t .t.) 


io’*if’ "o’ "f ' 


1 1 


o 




dt 


By requiring the trajectory to satisfy certain necessary 
conditions regarding path, equations (1.2.10) and (1.2.11), 
the Problem of Bolza has been reduced to the problem of 
minimizing J, a function of endpoints, subject to the i K 
algebraic constraints on the endpoints. 

Before proceeding with the minimization of J, it 
is appropriate to consider a graphical interpretation of 
the functional relationship for the state variables 
expressed in equation (2.1.9). Figure 2.1 shows a general 
state function y^(t ,y^ Q ,y^ f ,t Q , t f ) as a function of time. 
From the figure it can be seen that a change in the final 
state Ay^ while holding all of the other endpoints fixed 
causes a change in y^ for all values of t. Likewise a 
change in the final time At^, while holding all of the 
other endpoints fixed causes a change in the state y^ 
for all values of t. 

Formalizing this graphical interpretation in terms 
of differentials yields results which will be of value in 
the following sections. Using equation (2.1.9), the 
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differential of the state variables may be written as, 

dt f (2.1.22) 




dy. dy. dy, 

dt + ^tr dt o + 3t 


o 


f 


dy. dy 

+ dy_. Q dy jo + dy.. f dy jf 


Evaluating this expression at t = t^ gives 


^fVyjo’yjf’VV = dy 


ay± 


+ 


Sy ± 

dt 


o 


dy ± 

dt ° + 3t^ 


if dt 
dy 

dt _ + 
f dy 


dt 


I f 


f 


(2.1.23) 


l 


jo 


dy. 


dy jo + cFy 


jf 


dy 


jf 


•f 


The sum represented by the last term in the above equation 
can be separated into those products for which i ^ j and 
that for which i = j . Transposing dy^ to the right hand 
side, equation ( 2 . 1 . 23 ) becomes 


0 


3y ± 

dt 


dy. 


1 


f ^ 


dt „ + 


Sy i 

5y 


JO 


dy 


J'o 


dy. 


+ 


x 


dy 


jf 


dy if + 
f 

'dy ± 

_ i 

_ ay jf 

f 

i/j 




dy jf (2.1.24) 
i-j 


In the above equation the repeated subscripts on 
the last term do not imply summation. Since t, t Q , t^, 
y^ Q , and y^ have been assumed to be independent, equation 
(2.1.24) implies that ■ 
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<5t~ f 


a n 


By evaluating equation (2.1.22) at t = t and 
following arguments similar to the ones above, it can be 


shown that 


I = a n 

; o ^ 


(2.1.31) 




(2.1.32) 
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= o 




dy. 


oy 


jo 


o 

i=J 


1 


(2.1.33) 


(2.1.34) 


These identities will be useful in the proofs of 
necessary and sufficient conditions in the next section. 


2.2 Derivation of Transversality Conditions 

In determining the functional relationships in 
the last section* it was assumed that the control and 
adjoint variables were chosen so as to satisfy the Euler- 
Lagrange Equations (1.2.10) and (1.2.11). Equations 
(1.2.10) and (1.2.11) are referred to as the first path 
necessary conditions. In this section the endpoint neces- 
sary conditions (transversality conditions) are derived 
assuming that the first necessary conditions for path are 
satisfied . 

The solution to the path necessary conditions 
determines one or more trajectories (see section 1.3)* any 
of which may be expressed functionally as a set 


t ,t )] as shown in section 2.1. Once the functions 
representing one of these trajectories Is substituted into 


the Integral in equation (2.1.21), the Integration can be 
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performed. It is therefore clear that once the trajec- 
tory is specified , J is a function of only the parameters 
y^ Q , t Q and t^. Specifying the path reduces the 

problem of minimizing J to the well-known problem of 
finding the minimum of a function of several variables 
subject to algebraic equations of constraint (Bryson and 
Ho , 1969 ). 

It Is shown in Appendix A that if the arguments 

of J in equation (2.1.21) are to satisfy the constraints 

and minimize J, then it is necessary that the partial 

* 

derivatives of the auxiliary function, J shown below, 
with respect to y^ Q , y^ f , t Q , and t f all be equal to zero. 

The J function is defined by J = J + p^ijj ^ where 
J is given by equation (2.1.21). Using the definition of 
the function G from equation (1.2.8), J may be function- 
ally represented as 


U y if ^o^f 





ay i‘ 

wr _ 


dt 


( 2 . 2 . 1 ) 


In the above equation it Is understood that y^, and 

u. are all functions of the set ft,y. ,y.„,t ,t„). In 
writing the functional relationship shown above, It has 
been assumed that the controls u^. have been chosen in an 
optimal fashion In accordance with the control variable 
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Euler- Lagrange Equation (1.2.11). This is indicated hy 
the superscript o on u^ and on H. The partial deriva- 
tive of J with respect to y^ Q can now he written: 


lyT + 


i d y ± d 2 y ± 

r. 1 ~Sy~. dt + ^i c3y . ~dt" 
jo ^jo °Jo 


( 2 . 2 . 2 ) 

Here Leibnitz Rule ( Hildebrand > 1948 , p. 360) has been 
used for differentiation of an integral with respect to a 
parameter. Using the identity 


dt A i "dy. 


a y ± dx j _ 5y i 
^i dy . dt' ^ dt "SyT”" 

■jo J jo 


(2.2.3) 


and expanding -v- , equation (2.2.2) may be written as 

jo 


3h 8h 3 n 


(2.2.4) 


sh 3u k 3 n 3 n 3 h 3 n 

5y~3ir 


dt + X 


Terms under the integral sign may be combined to give 


dj _ du 

dy . ~Sy . ^i 

J 30 J jo 


- L 


(2.2.5) 


sh 3 m 3x i 


dx. dt 


-j - dH + WA i 
r .io ^r- 


, a H 3u k 


Sx, \ ay. 
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Note from equations (2.1.9) and (2.1.10) that once the 

optimal endpoints have been selected , 

By. dy. Bx. d\. 

dt~ a wr ~ dt • 


The integral term vanishes, since equations (1.2.3) * 
(1.2.10) and (1.2.11) were used to generate the func- 
tional relations (2.1.9) and (2.1.10). 

Using equations (2.1.27), (2.1.33) and (2.1.34) , 
it can he concluded that the sums represented hy the two 
remaining terms not containing G in equation (2.2.5) 


reduce to a single term,, - X 


With these considera- 


t 


o 


tions, equation (2.2.5) reduces to 


3J 


* 


^jo 


sg 

~5y 


O'o 


h 


= 0 


( 2 . 2 . 6 ) 


-x- 


By taking the derivative of J with respect to y^ and 
using arguments similar to those just presented (in 
this case equations (2.1.32), (2.1.28) and (2.1.29) must 


he taken into account), it can he shown that 


Bj 

"By 


jf 


bg 


+ \j 


0 


(2.2.7) 


Two more necessary conditions remain to be 
derived. These result from taking the partial deriva- 
tives of J with respect to the remaining two variables , 
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t Q and t^. Performing the first of these operations 
yields 


dj* __ dG 

’ST” dt~ 
o o 



dH 

?t o 


dX ± dy ± dy ± 

+ dt dt -- + x i dtot 
o o 


dt 


- H + X 


Sy i 

i dt 


( 2 . 2 . 8 ) 


o 

Here again Leibnitz Rule has been used; this time the 
limits of integration are functions of the differentiating 
variable. Using the identity 


d 

dt 


dy. 


X i dt 


1 


o J 


, dy i aX i Sy i 
- x i dtdt~ + dt~ dt~ 


( 2 . 2 . 9 ) 


o 


o 


and expanding ^jr-> equation (2.2.8) may be written as 

o 


d J* dG + 
^ t o 


-t , 


t 


dH 


dy. 


dH 3X i 


dy7 dt dT7 dt~ 

" x o x o 


o 


... .. s h 3 n 

du, eft d t ' d t d t "d t " 

k o o o 


dH dU k 


dt 


dy ± 
"i dt 

o 


t 


o 


dy . 

" H + X i dt 


x 


t 


o 


( 2 . 2 . 10 ) 


Terms outside the integral may be evaluated at 


the endpoints indicated and terms under the integral sign 
combined to give 
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dj* dG , TT 
3t“ “ + H 

o o 


dy 


t - x i 
o 


1 


t 


dt 


o 


t 


o 


+ X ± 


3y± 

- x i 


, dt 

t „ o 

, dt 

t o 

f 

tf 

o 


( 2 . 2 . 11 ) 


t 


o 


+ 



r 

( dH 

dy . \ 

X 

L 

- 

l 6y i 

+ 


o 


o 


dH 

ay± ) 

SX i 

I dH 

du, 1 
k 

di ± 

dt 

6t o 

du, 

1 kj 

^ _ 


The integral again vanishes identically for optimal paths. 
Using equations (2.1.26) and (2.1.30), the three terms 
outside the integral representing summations can also be 
equated to zero. With these observations, equation 
(2.2.11) reduces to 


dj* _ dG 

dt dt 
o o 


+ H 


t 


o 


= 0 


( 2 . 2 . 12 ) 


* 

By taking the derivative of J with respect to 
t f , following a line of reasoning similar to that just 
given, and using equations (2.1.31) and ( 2 . 1 . 25 ), it can 
be shown that 


, "5f 

dJ = dG 
'St t 



H 


0 


(2.2.13) 
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These results are summarized in the following 
statement ; 

2.2.1 Transversallty Necessary Condition for End - 
points . If a trajectory satisfies the Euler- Lagrange and 
state variable differential equations , equations (1.2.10)., 
(1.2.11) , and (1.2.3), and if the set E = fyich^if " ^o' 
t^,^] satisfies endpoint equations of constraint (1.2.4) 
and provides a local minimum of J with respect to small 
allowable variations in the endpoints, then the set E 
must satisfy equations (2.2.6), (2.2.7), (2.2.12) and 
(2.2.13) • 

These latter equations are referred to as the 
endpoint necessary conditions or, classically, as the 
transversallty necessary conditions. 

2,3 Derivation of Endpoint Sufficiency Conditions 

In the last section the function J was shown to 
be a function of the endpoint variables y,. o ,y ; . t Q , and 
t^ when evaluated along an optimal path. The function J 
is constrained, however, through the p equations of con- 
straint ijf ^ of equation (1.2.4). Sufficiency conditions 
for determining the minimum of a function whose arguments 
must satisfy algebraic equations of constraint are well 
known (Vincent, 1969 ), and, for reference, the sufficiency 
conditions are derived in matrix notation in Appendix A. 
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Before presenting a statement of the sufficiency 
condition, a brief discussion and definition of notation 
are in order. Since the algebraic equations of con- 
straint for the Problem of Bolza define relationships 
among the endpoint variables, the endpoint variables are 
not all independent. Since there are p equations of con- 
straint and (2n + 2) endpoint variables, there are only 
(2n - p + 2) independent endpoint variables. The p 
dependent variables are determined by the p equations of 
constraint. Any p of the variables can be considered to 
be the dependent variables. The choice is one of con- 
venience. Let the p dependent variables be denoted by 
the column vector w and the remaining (2n - p + 2) inde- 
pendent variables be denoted by the column vector _v. 

Let the vector represent a vector whose elements are 
the f constraint functions. Equation (2.3.1) summa- 
rizes these relations . 



t 1 (w,v)' 

. w = 




( 2 . 3 . 1 ) 


q = 2n - p + 2 

The identification of the elements of v_ and the 

elements of w with the endpoints y. ,y. _,t , and t„ is 
— J xo ^xf o’ f 
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arbitrary except that the set of § equations must con- 
tain every element of w and, in addition, every equa- 
tion must contain at least one element of w. It is con- 
venient to define an additional column vector r_, whose 
first elements are the dependent variables and last 
elements the independent variables: 


w 


1 


r = 


w 


P 


v- 


v. 


(2.3.2) 


With these vectors define 


dJ 

drdr 


dJ 


as a (2n +2) by 
* 


(2n + 2) matrix with elements a^^ = . Let the 

I 3 


matrix $ be defined by 


cH 

-1 ; 

Bijf 

_ <5w 

■ 

dv 


(2.3.3) 


where 


cTw 


an 


is a p by p matrix with elements a^. = 


and 


d'j/_ 


is a p by q matrix with elements a. . 


Bn 


ij dv . ‘ 


It is shown in Appendix A that the $ matrix is the linear 
transformation which transforms differential changes in 
the independent variables into differential changes in 
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the dependent variables . The $ matrix has p rows and q 
columns. Finally, define the (2n +2) by q partitioned 
matrix Q as 



(2.3.4) 


where I represents a q by q identity matrix. 

With these definitions the endpoint sufficiency 
condition may now be stated: 

2.3.1 Endpoint Sufficiency Condition . If E 
represents a set of endpoints and multipliers [y io -,y^ f , 
t t n ] which satisfy the transversality necessary 
condition for endpoints, then a sufficient condition for 
the set E to represent a local minimum of the function 
J with respect to small allowable variations in the end- 
points is that the quadratic form 



in the differentials dv must be positive definite when 
evaluated at the stationary point E. 

To implement this sufficiency test, it is neces- 


sary to evaluate the elements of the matrix § and the 
elements of the matrix . Evaluation of elements 

of $ represents no problem since the functional form of 
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the constraints is specified in the problem statement. 
However, the analytic evaluation of the second partial 
derivatives of J with respect to the endpoints is not so 
simple . 

* 

The second partial derivatives of J can be 
obtained by taking the partial derivatives of the trans- 
versality necessary conditions with respect to the end- 
points r. . 



( 2 . 3 . 6 ) 




(2-3.9) 


where in the above equations i = 1,2,..., 2n+2. The func- 
tional form of G as a function of the endpoints is speci- 
fied by the statement of the problem. However, the 
functions \ . and H are not known functions of the end- 

t) 

points until the state variable and Euler- Lagrange differ- 
ential equations have been integrated analytically. 

Since analytical integration is often difficult 
or impossible, it would be desirable to evaluate the 
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partial derivatives of and H with respect to the end- 
points in terms of functional forms specified in the 
statement of the problem. A complete set of relation- 
ships of this type were not found. Unless future inves- 
tigators establish such relationships , analytic applica- 
tion of the sufficiency condition requires an analytic 
solution of the state variable and Euler- Lagrange differ- 
ential equations. 

Some interesting relations of this type are easily 
obtained however. Each of the elements of the matrix 

is composed of a sum of a second partial 

derivative of G and a second term. The matrix can there- 
fore be expressed as the sum of two matrices. 





( 2 . 3 . 10 ) 


where the matrix A is determined from equations (2.3.6) - 

-X- 2 

(2.3.9). Since J and G are of class C by hypothesis, 

•X* . . 

both J and G must be symmetric about their major diago- 
nals. The obvious conclusion is that matrix A must also 
be symmetric . By equating symmetric elements of A, the 
following identities can be established: 






( 2 . 3 . 11 ) 



45 






aH o 

3y lf 




( 2 . 3 . 12 ) 

( 2 . 3 . 13 ) 





3t dt ~ 

o I 


( 2 . 3 . 14 ) 


In addition., the following relations can be established 
by considering the functional relationships exhibited in 
section 2.1. 


3H f 3H , - ax if 

dt f 'St' 1 if Tt' f " 


( 2 . 3 . 15 ) 



analytical integration of the state variable and Euler- 
Lagrange equations poses an interesting problem for future 
investigations . 



2.4 Relations to Classical Theory 
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Bolza (1961, pp. 102-103) gives an excellent 
summary of the various classical approaches to the devel- 
opment of necessary and sufficient conditions for vari- 
able endpoint problems. The classical problem in the 
calculus of variations is to minimize the integral of a 
function P = F(x,y,y') unconstrained by differential 
equations of constraint. Here x is the independent vari- 
able , y is the independent variable, and y' represents 

dy/dx . The first and second order variation of the 

2 

integral are written as &J and 6 J, respectively, while 
<5y and 8x represent variations in the dependent and inde- 
pendent variables. 

Because of the pertinence of Bolza’ s remarks to 
this presentation, his historical synopsis is quoted in 
detail : 

Three essentially different methods have been 
proposed for the discussion of problems with vari- 
able end-points: 

1 . The method of the Calculus of Variations 
proper : It consists in computing 6J and <$2j 

either by means, of Taylor's formula or by the 
method of differentiation .with respect to e , . . . 

and discussing the conditions <5J ■ = 0, <$2j > 0 . 

The method was first used by LAGRANGE . . . [. (l867» 
pp. 338, 345)3. He gives the general expression for 
6 J when the endpoints are variable, viz.: 
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8 J 



&y 


F. 


y 


~~ F r 

dx y ' 


dx 


f 

+ [F6X + F r 6y] 

^ o 


(2.4.1) 


and derives the conditions arising from 6J = 0. 

The second variation for the case of vari- 
able end-points was first developed by Erdmann 
. . . [(1878, p. 364)]. He finds 


-x 


§ 2 J 


f R(u6y' - u r 5y) 2 dx 


(2.4.2) 


x. 


u 


o 


rh 


F8 2 x + F yl 6 2 y + 2F y , 6x8y + 2F y ,8x5y' 


dF c 2 


u 1 


+ ^ 6X ^ + p y'y +F yVF 6y 


0 


where u is an integral of Jacobi’s differential 
equation . . . [ (Bolza, 1961, p. 49) ] • By con- 
sidering such special variations for which 6y = 
CUj he makes the integral vanish and thus 
reduces the question to the discussion of the 
sign of the remaining function of the variations 
6x.j 6y., 5 2 x. , 6 2 y. . These variations are con- 
nected Toy relations 1 which depend upon the spe- 
cial nature of the initial conditions .... 

For the general integral 

x f 

J = / F(x^y 1 ,yg. . . . ,y n Jy 1 < ,J 2 ' , ■ ■ ■ ,r n ' )dx 

x o (2.4.3) 

where y-, . . . ,y are connected by a number of 

finite or differential relations, the second 
variation in the case of variable endpoints was 
studied by A. Mayer . . . ((1896, p. 436) 3 ■> for 
the integral in parameter-representation 
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f 

J = J F(x,y,x f ,y ' )dt (2.4.4) 


by Bliss . . . [(1902, p. 132)]. 

2 . The method of Differential Calculus : 

This method is explained, in a general way by 
Dienger . . . [(1867)]. It decomposes the prob- 
lem into two problems by first considering 
variations which leave the end-points fixed, 
and then variations which vary the end-points, 
the neighboring curves considered being them- 
selves extremals. The second part of the prob- 
lem reduces to a problem of the theory of 
ordinary maxima and minima. This method has 
been used by A. Mayer in an earlier paper on 
the second variation in the case of variable 
end-points for the general type of integrals 
mentioned above . . . [(Mayer, 1884, p. 99)]- It 
is superior to the first method not only on' 
account of its greater simplicity and its more 
elementary character, but because-~by utilizing 
the well-known sufficient conditions for ordi- 
nary maxima and minima- -it leads, in a certain 
sense, to sufficient conditions If combined 
with Weierstrass 1 s sufficient conditions for 
the case of fixed end-points .... 

3- Kneser's method : This method, which 

has been developed by Kneser . . . [(1900)], is 
based upon an extension of certain well-known 
theorems on geodesics. It leads in the simplest 
way to sufficient conditions, but must be sup- 
plemented by one of the two preceding methods 
for an exhaustive treatment of the necessary 
conditions .... 


Later investigators of variable endpoint prob- 
lems, e.g., (Householder, 1937), (Bliss, 1946), and 
(Hestenes, 1966) have followed the first method quoted 
from Bolza, "the method of the Calculus of Variations 
proper." To the best knowledge of the author, there 
have been no further developments or exposes using 
"the method of Differential Calculus" since that of 
Bolza (1904) . 
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The method developed in this presentation is 
essentially the second method, "the method of Differen- 
tial Calculus," quoted above from Bolza. Paticular 
attention should be paid to his remarks concerning this 
method. These remarks are consistent with fundamental 
propositions of the last three sections, namely: 

A comprehensive sufficiency condition for 
variational problems with variable endpoints 
is obtained by applying two independent tests. 

Test A and Test B below, each of which is 
applied separately. 

Test A. Satisfaction of the endpoint suffici- 
ency condition given that the trajectory 
satisfies necessary path conditions. 

Test B. Satisfaction of the path sufficiency 
conditions for the endpoint fixed. 

In addition, Bolza 's remarks indicate that 

For a solution to be optimal, it is necessary 
and sufficient that its endpoints satisfy the 
endpoint sufficiency condition (Test A above) 
and that its path satisfy the fourth necessary 
condition with its endpoints considered fixed 
(Test B above ) . 

Although this presentation has not undertaken a rigorous 
proof os this hypothesis, the examples and analysis have 
given every indication that the hypothesis is valid. 


2.5 Geodetic Example 

As an example of the application of the suffi- 
ciency condition for endpoints, consider the problem of 
determining the minimum distance from the origin to any 
point on a parabola of the form 



2 , T 

y « x + b 


50 

(2.5.1) 


In control notation the problem may be formulated as 
follows : 

Minimize 

f 

J = J ds (2.5.2) 

s 

o 

subject to the state variable differential constraints. 

If = cos g, (2.5.3) 

§J = sin g, {2 - 5 - J +) 

and endpoint constraints , 


^0 = 


x = 


s = 


0 , 

0 , 

0 , 


y f -x f -b 


= 0, 


(2-5-5) 

(2-5-6) 

(2-5-7) 

(2.5.8) 


The angle g is the angle between the positive x axis and 
a tangent to the curve . Here x and y are the state vari- 
ables , g is the control variable, and s is the independent 
variable analogous to t in the formulation of earlier 
sections . 

2,5-1 Necessary Path Conditions . The H and G 
functions are 

H - X cos g -f X sin g - 1 

x y 


(2-5.9) 




2.5-2 Functional Relations . Integrating the 
state variable equations (2.5-3) and (2. 5-^-) with the opti- 
mal constant control g between the general Initial point 
(x o ,y o ,So) and general final point (x^y^jS^) results in 

X f - X Q = (s f - s Q ) cos g (2.5.17) 

y f “ y© = ( s f - S Q ) sin s 


(2.5.18) 
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Solving for the control 

y r - y n 

tan g - -S (2.5-19) 

x f “ x o 

Squaring both sides of equations (2.5-17) and (2.5.18) 
and adding yields the identity 


( s f " S Q ) 2 = (x f - x q ) 2 + (y f - y Q ) 2 (2.5-20) 

Solving equations (2.5.17) and (2.5.18) for the controls 
gives 

cos g = -i (2.5.21) 

s f “ s o 

and 

7-p " 

sin g = — — y— ~ (2.5.22) 

S f S o 

Since the control is constant, the control is not a func- 
tion of the independent variable in this case. For other 
problems the control may be a function of the independent 
variable as well as the endpoints. 

Integrating the state variable equations again 
between the general initial point ( x 0 >y o > s 0 ) and general 
intermediate point (x,y,s) and substituting the optimal 
control from equations (2.5.21) and (2.5.22), and 
rearranging yields 

x “ x o + vr^r ( s - s o> (2.5.23) 

f o 



(2.5.24) 
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It is seen from the above equations that the state vari- 
ables are clearly functions of coordinates of the initial 
and final state variables and of the initial and final 
values of the dependent variable. 

The first integral of the Euler- Lagrange equa- 
tions is 

X cos g + X sin g - 1 = 0 (2.5.25) 

a y 

Solving this equation with equation (2.5.13) for X x and 
X and observing equations (2. 5. 21) and (2.5.22) gives 


= x f x o 
s f - s o 


x f - x o 


(x f -x o ) 2 + (y.p-y n ) 2 


(2.5.26) 


'f “ O ' 


y. 


V 


\ = 


o 


f 


- s 


o 



(2.5.27) 


Two forms are given above for the Lagrange multipliers 
as functions of endpoints; either is correct. If the 
second set is used, the J function will be independent 
of s^ and s Q . In either case it is clear that the 
Lagrange multipliers can be written as explicit functions 
of the coordinates of the initial and final states and of 
the initial and final values of the dependent variables . 
Equations (2.5.19), (2.5.23), (2.5.24), (2.5-2 6 ), and 
(2.5.27) bear out the functional dependencies hypothesized 
for control, state, and adjoint variables in section 2.1. 



54 


Note that in deriving these equations , only path neces- 
sary conditions have been used. The transversality 
necessary conditions for endpoints have not been -used. 


2 . 5.3 Necessary Endpoint Conditions . The trans- 
versality conditions, (2.2.6), (2,2.7)* (2.1.12), and 


(2.2.13) yield the following equations 
* 


5j 

5 y~ ti 2 


1 


8j 

5 x 

8J 

cTs 


yo 


^3 " x xo 


0 


0 


o 

* 


M4 


+ H o = 0 


5 j 

= li-L + X 

s * 

5 j 


yf 


= 0 


55 c] 

5 j 
5 s] 


- 2p 1 x f + x 


xf 


- 0 


* 


(2.5.28) 

( 2 . 5 . 29 ) 

( 2 . 5 . 30 ) 

( 2 . 5 . 31 ) 

(2.5.32) 


= - X 


xf 


cos g f - x^ r .p Sin g. p + 1 = 0 (2.5.33) 


yf 


Since the initial point is fixed, the initial point trans- 
versality equations ■ give no useful information. 

To find the optimal endpoints, eliminate 
between equations (2.5.31) and (2.5-32), yielding 

X xf -t- 2 X yf x f = 0 ( 2 . 5 . 34 ) 


Substituting X x ^ and X ^ from equations ( 2 . 5 - 26 ) and 
(2.5.27) into equation ( 2 . 5 . 34 ) yields 
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(2.5*35) 


Finally, multiplying through by the radical and imposing 
endpoint constraints ( 2 . 5 - 5 ) and ( 2 . 5 - 6 ) gives 

x f (l + 2 y f ) = 0 ( 2 - 5 - 36 ) 

The necessary conditions are satisfied if either term in 
the above equation is equal to zero. Solving equation 
( 2 . 5 . 36 ) and equation ( 2 . 5 . 8 ) simultaneously gives the 
two solutions 

x f = + ~ - b y-f = - (solution A) 

and > ( 2 . 5 - 37 ) 

x„ = 0 y„ = b (solution B) 

1 1 ) 

These endpoints and the corresponding multiple 
solutions for b < - 75 are shown in Figure 1-3 on page 14. 
From the symmetry of the parabola, it is expected that 
either the plus or the minus sign in equation ( 2 . 5 * 37 ) 
will determine a solution giving the same value of dis- 
tance. For this reason a distinction has not been made 
between the two. The necessary conditions used so far 
have provided no means for determining under what cir- 
cumstances solution A (or solution B) is the optimum. 

In this case of multiple stationary solutions, the .endpoint 
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sufficiency condition will provide a means for deter- 
mining the true optimum. 

Before examining the sufficiency conditions , the 
parameter will he evaluated in terms of the general 
endpoints for future reference. From equations (2.5-27) 
and (2.5.31) it is observed that 




1 


1 


yf 


y f-y 0 

(Xf,-X n ) 2 + (y.p-yj 2 


(2.5.38) 


"f "o' ' Wf "O' 

2.5.4 Sufficiency Endpoint Condition . To evalu- 
ate the endpoint sufficiency conditions, it is instruc- 
tive to first determine the $ and Q matrices of equations 
(2.3.3) and (2. 3 .4). The constraints are 


*1; 

I' 2 ; 

*3: 


7r 


0 


x = 0 
o 


y f - x f 


o 


b 


D 


0 


(2.5.39) 

(2.5.40) 

(2.5.41) 

(2.5-42) 


Since there are four equations and six endpoints, there 
are two degrees of freedom. For conveniency let and 
s^ be the independent variables and y o ,x Q ,s o , and y^ 
be the dependent variables. Then in the notation of 
section 2.3 
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v = 



y o 


y o 


y o 

x^ 


X 


X 

o 


0 


o 

s 

. r = 

s 

, i = 

s 

o 


o 


o 

y f 


y f 


2 , 

y f - Xf -d 


x f 

«* 


s f 



(2.5.43) 


Evaluating the matrix of partial derivatives of with 


respect to the independent variables gives 



0 0 

0 0 

0 0 

-2x f 0 


(2.5.44) 


Evaluating the matrix of partial derivatives of jr_ with 
respect to the dependent variables gives just the iden- 
tity matrix 


c5w 


10 0 0 
0 10 0 
0 0 10 
0 0 0 1 


(2.5.45) 


The inverse of this matrix is obviously the identity 
matrix. From equations (2.5,44) and (2.5-45) the $ 
matrix can be computed 
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■ 81 "j 

-1 

" 

dw j 


dv 


o o 

o o 

0 0 

2x f 0 


(2.5.46) 


The Q matrix is formed by adjoining to $ an identity 
matrix with the dimensions equal to the number of inde- 
pendent variables. In this case there are two indepen- 
dent variables . The Q matrix is 


Q = 


0 

0 

0 

2x i 

1 

0 


0 

0 

0 

0 

0 

1 


. (2.5.47) 


With the 


•fr -1 


dj_ 
dr 6 r 


not yet evaluated, the endpoint suffi- 


ciency condition reduces to the condition that 


[dx f ds f ] 


4x 


4x^ 


-O * 

2 a^j 

"a 2 
oyf. 

3 2 j* 

ay f ax f 

a g j* 

dx 2 


“I** 


2x 


"f* 


f ds f 6y f 

d 2 J* 

ds f dx f 


must be positive definite. 


2x 


a 2 / 

f 

■x2 * 

afj 

ax f a s f 


~b 


a 2 j* 
8s 2 


dx. 


ds 


f 


(2.5.48) 
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If J* can be written so that it Is not a function 
of s^, the sufficiency condition will he reduced to a 
simple inequality involving dx^ only. From the trans- 
versality equations (2.5-31) - (2.5-33) and the func- 
tional relations for X x and A. , equations (2.5.2 6) and 
(2.5.27 ), it is seen that this can be done. 

Therefore, the sufficiency condition reduces to 
the condition that 

+ <>° 


2 . 

Since dx^, is always positive, the expression 111 paren- 
thesis must be positive in order to satisfy the suffi- 
ciency condition: 


a 2 j* + 4x ^ 2 j!_ + aV > 0 
a y f f axl 


This result is identical to the result that would 


have been obtained if the fixed endpoint coordinates, y , 

x , and s had been excluded from the G function. This 
o o 

situation is similar to the transversality necessary con- 
ditions in that the initial points yield no information. 
From this example and previous experience with endpoint 
conditions, the following conclusion is drawn: Wo useful 

information concerning either necessary or sufficient 
conditions results from Including in G constraints which 
merely fix a given endpoint coordinate. 
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Substituting the functional forms for X ^ 


and 


'X ^ not involving s^ from equations (2.5.26) and ( 2 . 5 - 27 ) 

* 

into the first partial derivatives of J with respect to 
y f and x f in equations (2.5.31) and (2.5.32) gives 


c)j' 

6y f U 1 


y f - y Q 


(x f -x o ) 2 + {y r -y n ) 2 


(2.5.51) 


f “ O ' 


6j 

cFxr 


2^ l x f + 


x f~ x o 


(v-v ) 2 + (v-yn ) 2 


( 2 - 5 - 52 ) 


if ~o' ' wf j o - 

Forming the required partia-l derivatives results- in 


^ 2 -x- 

6 J q , ^f 

5 = - 2|il + 

OX^. 


( x c~x n ) 2 + (y f -y 0 ) 2 " x -p( x r" x o) 




D 


( 2 . 5 . 53 ) 


•n2 * 

5 J 

^ y f^ x f 


(x f -x 0 )(y f -y 0 ) 

D 


ah* ( x f - x 0 ) 2 + h f -y 0 ) s -y f (y f -y 0 ) 

% 2 D 


( 2 . 5 - 54 ) 


(2.5.55) 


^2 * 

6 J 


(y f -y 0 ) ( x f - x 0 ) 

D 


( 2 . 5 . 56 ) 


where 


2 1 


d = [ (x f -x D ) + (y f ~y 0 ) ] 2 


(2.5.57) 


Comparing equations (2.5.54) and (2.5.56) verifies the 


symmetry of the 


S 2 f 

drdr 


matrix. 



6l 


Evaluating these derivatives using initial point 
constraint equations (2,5.5) and (2.5.6) and from 
equation (2.5-38) and substituting them into the endpoint 
sufficiency condition, equation (2.5-50) gives 

4 

4 x„ 


4 x 


■f 




Vf 


2 " 2x3/2 
(x f + y f ) 


+ 2 




x f + ^f 


(2.5.58) 




y f 

(x f 2 + y f 2 ) 3/2 


> 0 


For solution B (x^ = 0 , y^, = b) this condition reduces to 


^ + 2 > 0 
b 


( 2 - 5 - 59 ) 


From the geometry in Figure 1 . 3 * it can be seen that 
b is negative.; The" -condition therefore requires that 

0 > b > - (2.5.60) 

Then solution B as . shown in Figure 1.3 is optimum. 

r f 


For solution A .(x.p = + -\/-~~b , y. 


1 


‘■f - - V 2 ~ ' •'f ' 2 

the end sufficiency condition (2.5.58) becomes 

4 b 2 + 2b + 


), 


(- i - *) 3/2 


(2.5.61) 


VT 


T 

% 


> 0 
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Combining terms yields 


4b 2 + 3b 

(- i - 6 



> 0 . 


(2.5.62) 


In order for the denominator to be real 

b < - ( 2 . 5 . 63 ) 

Under this condition inequality (2.5.62) is satisfied 
only if 

b < - | . (2.5-64) 

Therefore solution A shown in Figure 1.3 is optimum for 
b less than - 75. The optimal solution is summarized below. 

x f « 0, y f = b 0>b>--| (2.5.65) 

x f = + y-|-b, y f = - b < - ~ (2.5.66) 


This simple example has been analyzed in great 
detail to emphasize the concepts developed in earlier 
sections and to reinforce and illustrate the notation. 

A more complex example of the endpoint sufficiency condi- 
tion is presented in Chapter 5- 
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2.6 A Numerical Algorithm 

In order to apply the endpoint sufficiency condi- 

* 

tion, the matrix of second, partial derivatives of J with 
respect to its arguments must he determined. From equa- 
tions (2.3.6) - (2.3.9) it is seen that each of these 
second partial derivatives is composed of two terms. The 
first term., in all cases, is a second partial derivative 
of the function G. This derivative can be computed 
analytically from information given in the statement of 

the problem. The second term of each second partial 

*)(■ 

derivative of J can be written in one of the following 
forms : 

( 2 . 6 . 1 ) 
( 2 . 6 . 2 ) 


dM 


t=t J 




dM 


t=t 


o 


31 


if 


or 




dM 


t-t. 


o 


dr. 

10 


(2.6.3) 


(2.6.4) 


where M represents any of the quantities H, X-j, Xg, . 
X n and 


• • > 
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r represents any of the state variables y^ or the 
independent variable t . 

These derivatives cannot be evaluated analytically with- 
out obtaining an analytic solution to the set of state 
variable and Euler- Lagrange differential equations. For 
most problems of practical interest in the calculus of 
variations,, the set of nonlinear state variable and 
Euler- Lagrange differential equations cannot be inte- 
grated analytically. Therefore , the implementation of 
the endpoint sufficiency condition in most cases requires 
the numerical computation of partial derivatives of the 
forms expressed in equations (2.6.1) - (2.6.4). 

Fortunately, this is not conceptually difficult 
for most problems in engineering which have separated 
end constraints. End constraints are separated if none 
of the endpoint constraints involves both initial values 
and final values; the constraints always relate initial 
values to other initial values, or final values to other 
final values . 

The function M evaluated at t = t will be indi- 
cated by a subscript o: 

M o = ’W (2.6.5) 

The function M evaluated at t = t f will be indicated by 
a subscript f: 



65 


M f = M f ( y io ,y if ) ' 

Before the sufficiency condition test is applied, 

the problem is first solved using the necessary condi- 

*)(- 

tions yielding nominal endpoints y^ Q , t Q , anc ^ 

-X" "X - 

and nominal Lagrange multipliers and F°r brev- 

* -x- * 

ity, let _r Q represent a vector with elements (y^j y 2 0 > 

* * . * _ 

. . . , y nQj t Q ), and represent a vector with elements 

4c- -X- -x- -X- -X- 

t , y"2f 5 • • • j y-nf’ and M be a function evaluated 

with the nominal endpoints . 

Numerically the derivative (2.6.1) can be approx- 
imated as 


cSM f 

cir\. 

10 


«f(Uo’ 


* 

"20 ’ 


> r io 


-I- A, 


> Lf) 


* 


M, 


A 


(2.6.7) 


where A is a small change in the nominal initial vari- 
able r^ Q . If the state variable and Euler- Lagrange equa- 
tions are then numerically integrated forward with the 
nominal Lagrange multipliers , the final nominal endpoint 
will not be reached. The n initial Lagrange multipliers 
must be adjusted in order to obtain the final nominal 
endpoint again. Since the n initial Lagrange multi- 
pliers give only n degrees of freedom, the nominal end- 
point can be reached only if M is a function of n or 
less than n Independent final values. This will be true 
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if there is at least one equation of constraint involving 
the final values. With these new multipliers, the differ- 
ential equations are integrated forward to the final point 

-X- 

is then evaluated from the resulting final Lagrange 
* 

multipliers and r^.. With evaluated, the desired par- 
tial derivative can be evaluated using equation (2.6.7) • 

The derivative (2.6.2) can be approximated 
numerically as 


SM 0 


M, 


( *. 
~o -o * 


* 

'if • 


* 

? 2f J 
A 


* 

' r if + A > 


. . ) - M, 


* 


o 


( 2 . 6 . 8 ) 


In the above equation, M q is evaluated by making a small 
* 

change in r^„, while leaving all the other values 
unchanged. A set of final Lagrange multipliers is then 
determined so that a backward numerical integration in 

■3fr 

time will yield the nominal initial values r . The 
J —o 

quantity M q is evaluated using the resulting initial 
Lagrange multipliers and r Q . With M Q computed in this 
manner, the desired partial derivative can be evaluated 
using equation (2.6.8). 

The derivative (2.7-3) can be approximated numer- 
ically as 




or 


f ^ r lf’ r 2f’ r if + Aj M f (2.6.9) 

if A 


Here, is evaluated by making a small change in the 
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* 

nominal final point coordinate r^, while leaving all of 

■* 

the other final coordinates and the initial point r^ 

unchanged. A set of initial Lagrange multipliers is then 

determined so that a forward integration from the nominal 

. * 

initial point will yield the varied final point (r^; 

? 2 f’ • • • > r if + A; ...). The forward integration is then 

performed to the varied final point * and is evaluated 
using the resulting final Lagrange multipliers and the 
coordinates of the varied final point. 

The final derivative (2.6.4) can he approximated 
numerically as 



; r . + A; . 

xo ’ 


( 2 . 6 . 10 ) 


Here> M Q is evaluated by making a small change in the 

nominal initial point while leaving all of the other 

initial coordinates and the final point unchanged. A 

set of final Lagrange multipliers is then determined so 

that a backward integration in time from the nominal final 

point will yield the varied initial point (r^ Q ; r 2 0 > 

■x* x 

r^ Q + A, . . . ) . The backward integration is then performed 
to the varied initial point; and M q is evaluated using 
the resulting initial Lagrange multipliers and the coordi- 
nates of the varied initial point. 

Using the above techniques; the matrix of second 

■X' 

partial derivatives of J with respect to its arguments 



68 


can be evaluated. Because of the identities (2.3-11) - 
(2.3-14), there is some choice as to which of the above 
derivatives is used to evaluate the sufficiency condi- 
tion. It is a simple matter to numerically evaluate the 

matrix Q from the nominal initial and final points and 

2 

to test the matrix ^ for positive- def initness . 

The details of programs written to perform these opera- 
tions are left for discussion in Chapter 3- 

Selecting the correct Initial Lagrange multi- 

/ 

pliers so that the desired final points are reached as a 
result of integration is termed a problem with mixed end 
conditions, or a two-point boundary value problem. The 
numerical implementation of the sufficiency condition 
depends strongly upon the existence of numerical tech- 
niques for solving two point boundary value problems. 
These techniques are explained in detail in the next 


chapter . 



CHAPTER 3 


NUMERICAL SOLUTION OF TWO-POINT BOUNDARY VALUE PROBLEM 

Whether one attacks optimal control problems from 
the point of view of the classical calculus of variations 
or by application of the Maximum Principle, one invariably 
must solve a set of ordinary differential equations sub- 
ject to both initial and final boundary conditions. Such 
problems are referred to as two-point boundary value 
problems. The theory of ordinary differential equations 
subject to initial conditions has been well developed, 
both analytically and computationally, for some time. 
However, the theory of two-point boundary value problems 
has been slower to develop, especially from a computa- 
tional point of view. In the past deca.de, however, signifi- 
cant advances have been made in the application of large 
scale digital computers to the solution of two-point 
boundary value problems. The primary motivation for these 
advances has come from the study of optimal control 
theory (Handelsman, 1966; McGill and Kenneth, 1964) . 

This chapter discusses the use of a numerical 
technique for solving two-point boundary value problems 
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called the generalized Newton-Raphson method. Section 
3-1 briefly describes the relationship of this method to 
other numerical algorithms currently available. The 
method is explained in section 3*2 and a brief discussion 
of the actual computational procedure is given in section 
3.3. 

3.1 Numerical Methods Available 

Currently,, there are many numerical algorithms 
for solving two-point boundary value problems; the cur- 
rent question is one of deciding which to implement. To 
aid in this decision a. brief comparison of available 
computational algorithms is now made . The point of view 
of Kalman (1964) is taken in the discussion which fol- 
lows . 

The computational solution of optimal control 
problems always requires that the following five condi-? 
tions be satisfied: 

(1) state-variable differential equations (1.2.3) 

(2) algebraic equations constraining state variable 
endpoints (1.2.4) 

(3) transversality algebraic equations (1.2.12) - 

(1.2.15) 

(4) optimal control conditions: max u H(y^,x^,u^, t) 
(equation (1.2.11) or (1,6.2)) 
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(5) adjoint variable Euler -Lagrange differential 
equations (1.2.19) 

The conditions above define a two-point boundary value 
problem. Five computational algorithms which attempt to 
satisfy these conditions by successive convergent approx- 
imations are discussed below. 

3.1.1 Flooding Technique . A great deal can be 
learned about the optimal control of a system without 
attempting to obtain a solution to the two-point boundary 
value problem. If there are p unknown initial values, 
the problem is first relaxed by not requiring p conditions 
of (2), which constrain a final state coordinate, to be 
satisfied. The problem then becomes an initial value 
problem with one or more arbitrary initial values. If 
the number of arbitrary initial values is one or at most 
two, the method of flooding (Vincent and Brusch, 1966) 
is practical. 

Solutions are generated satisfying conditions (1), 
(3), (4), and (5) and Initial conditions of (2). The 
equations of conditions (3) and (4) are combined, elim- 
inating the p parameters, to yield a cutoff function. 

Xj 

Arbitrary values for the unspecified initial variables 
are chosen, and the resulting initial value problem is 
integrated numerically until one of the zeros of the cut- 
off function is encountered, indicating the transversality 
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conditions have been met. Using this method, the state 
variable endpoints corresponding to the arbitrary ini- 
tial values simply fall where they may. A family of 
solutions is generated using the above technique by sys- 
tematically varying the initial parameters throughout 
their range. The family specifies the region of end- 
point space throughout which solutions are possible. 

In addition, a mapping is obtained between unknown ini- 
tial values and final endpoints . 

Since only a single integration is needed for 
each solution, this method is economical for parametric 
studies when compared with more sophisticated techniques 
described later. The technique economically generates 
such a large number of results that solutions to all 
points in endpoint space may be easily visualized and 
areas of unusual interest quickly discovered. From the 
resultant manifold of solutions, the solution to a par- 
ticular two-point boundary value problem can be readily 
approximated . 

The disadvantages of this technique are the 
following : 

(1) Exact solutions to particular endpoints are 

not obtained. It is desirable, in the case of 
multiple solutions, to compare directly two 
differing solutions to the same endpoint. 
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(2) Often no optimal endpoints exist for certain 

ranges of initial values of the Lagrange multi- 
pliers (Vincent and Brusch, 1966 , pp. 27-28). 
These ranges are difficult to predict and inte- 
grations with initial values within these ranges 
yield no useful information. 

3.1.2 The Classical Indirect Method . To use 
this procedure , a series of solutions satisfying condi- 
tions (1), (4) and (5) are generated (see pages 70-71). 
This series converges to a solution satisfying conditions 
(2) and (3) • The usual method for satisfying conditions 
(2) and (3) is to numerically generate a matrix of partial 
derivatives of the final coordinates with respect to 
unknown initial values . This is done by making a sma.ll 
change in one of the unknown Initial values and observing 
the resulting changes in the final coordinates. By con- 
sidering the final coordinates to be linear functions of 
the unknown initial values, a simple matrix inversion 
yields new initial values and the iteration is repeated. 

The method requires good guesses for the initial 
values of the adjoint variables, particularly for high 
dimensional state spaces. The variations In the termi- 
nal boundary conditions with variations in the unknown 
initial values are often so sensitive that numerical 
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solutions are impractical. For further discussion, refer 
to Kopp and Moyer (1966, pp . 105-114) . 

3-1.3 Method of Gradients . This procedure is 
initiated by making an initial estimate of the control as 
a function of time; the corresponding trajectory is com- 
puted, and the value of the performance index calculated. 
The "direction" in state space for which the rate of 
change of the performance index is greatest is then deter- 
mined as a function of time. A small step is then taken 
in the gradient direction. That is, at each point in 
time, each state space coordinate is modified by a small 
quantity proportional to the "projection" of the gradient 
vector on its coordinate axis. The performance index 
is reevaluated and the iteration proceeds until the 
value of the index is stationary. As opposed to the 
indirect methods previously discussed which solve the 
first order necessary conditions on p. JO, this method 
is considered a direct method, since a direct search for 
the extreme value of the performance index is made. 

Methods previously discussed have been indirect since 
solutions have been generated using necessary conditions 
of the Calculus of Variations . 

In terms of the necessary conditions, listed on 
p. 70, conditions (1) and (5) are satisfied by each solu- 
tion of a series of solutions, and iteration proceeds 
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until (2), (3) and (4) are satisfied within some toler- 
ance. Convergence of this method does not depend upon 
good initial guesses for the controls. Although the 
method is relatively easy to program, its efficiency is 
limit ed, since the convergence slows in the neighborhood 
of the optimal solution. For a complete discussion of 
this theory, consult Kelly (1962, pp. 205-254). A 
detailed discussion of a numerical algorithm with an 
example is found in Hillsley and Robbins (1964, pp. 107- 
134) . 

3.1.4 The Second Variation Method . The second 
variation method (Kelly, Kopp, and Moyer, 1964) is also 
termed a direct method, since a direct search is made' 
on the functional. Again a series of solutions satisfy- 
ing conditions (l) and (5) is generated and iteration of 
the solution proceeds until conditions (2), (3), and (4) 
are satisfied. The main advantage of this method over 
the method of gradients is that convergence in the 
neighborhood of the solution is much improved. The sec- 
ond variation method determines the step size to be taken, 
while the gradient method provides no indication of step 
size. Again, convergence to a local minimum is independent 
of the initial solution estimates. Unfortunately, the 
computer time saved by this method is offset by computer 



76 

programming which is significantly more complicated than 
for the gradient method. 

3.1.5 The Generalized Newton- Kaphs on Method . 

This method^ which will he discussed in detail in the 
next section; is an indirect .method . It differs from 
classical indirect methods only in the technique used to 
solve the two-point boundary value problem. Necessary 
conditions (2) > (3); and (4) are satisfied by each solu- 
tion of a series of solutions; and an iteration is per- 
formed in function space which converges to a set of 
functions satisfying conditions (1) and (5) • The deter- 
mination of the step size is inherent in the method; and 
convergence is rapid near the optimal solution. Pro- 
gramming for this method is more complex than for the 
gradient method; but less involved than that required for 
the second variation method. 

The convergence of this method depends upon the 
choice of initial solution estimates. A sufficiency 
theorem for convergence has been developed by McGill and 
Kenneth (1963) for a general system of n second order 
ordinary nonlinear differential equations. Unfortunately; 
as the authors point out; problems of engineering interest 
seldom satisfy the stringent conditions of the theorem. 

The conditions of the theorem are sufficient; but not 
necessary; and many investigators; including this author; 
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have achieved convergence from elementary initial esti- 
mates which do not satisfy the conditions of the theorem 
(Moyer and Pinkham, 1964, pp . 91-106; Lewa.'llen, 1967). 
Sophisticated and effective algorithms have been devel- 
oped by Lewallen, Tapi ey, and Williams (1968) for 
stabilizing convergence of the iteration at some sacri- 
fice in efficiency. With these techniques , convergence 
envelopes have been achieved which permit initial esti- 
mates , generated with some of the initial values in error 
by over 100 per cent, to converge to the optimal solu- 
tion. If initial estimates are in too great of error, 
the iteration diverges. 

3.2 A Generalized Method of Newton- Raphson 

The Newton- Raphson method was first suggested by 
Hestenes (1949) for obtaining solutions to fixed end- 
point problems in the calculus of variations . The method 
was expanded and developed by Bellman and Kalba (1965) to 
include problems with a variety of boundary conditions. 
Bellman and Kalba referred to the method as quasilinear- 
ization. The method was applied to n-dimensional opti- 
mization problems formulated in state variable, control 
variable notation by McGill and Kenneth (1964) . The 
following discussion follows their work, except for the 
section on variable endpoint problems with final time 


free . 
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3.2.1 An Iteration for Nonlinear Initial Value 
Problems . Consider the following system of 2n nonlinear 
differential equations : 

Z = Z(Z;t) (3.2.1) 

subject to fixed initial point constraints 

Zj (0) = yj 0 J' = 1,2, ...,p (3-2.2) 

and final point constraints 

y k (t f ) = y kf k = 1,2, ...,2n-p (3.2.3) 

where 

F = (f!,f 2 ,...,f 2 n) (3.2.4) 

f i = f ‘ ‘ ,Y 2n ,t ) 1 = (3.2.5) 

A = (f i4 2 ^ • * * 4 2n _ p ) (3*2.6) 

If F is a function of t, P can be transformed into 
an autonomous system by observing the change <of variable 


n n +i = t 

(3*2.7) 

4n+l = 1 

(3.2.8) 

and adding the constraint 


y ( 2n+l ) o ~ 0 

(3.2.9) 


Thus, F may be considered as a function of only 
y without any loss in generality; such a functional rela- 
tionship is assumed in the following development. 
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Consider the truncated Taylor Series expansion 
of F(y) about some nominal state function ;y°(t) which 
satisfies initial condition constraints but does not 
necessarily satisfy equation (3.2.1): 


p lPc)] s p n°(t)i + 


OF' 


[y 1 (t) - z u (t)] (3.2.10) 


o 


o 


where 


z( t ) = [y 1 ( t ) 5 y ? ( t ) > • • • >y Pn (t)] 


2n' 


( 3 . 2 . 11 ) 


and 


r 0 F 

W 


is recognized as the Jacobian of the vector func- 
tion. Here y°(t) may be identified as an estimate to the 
solution of equation (3-2.1). If, at any time t, a new 
trajectory y^(t) is not equal to the estimate y°(t), 
equation (3-2.10) gives an approximate formula for esti- 
mating values of the vector function F [y^ (t ) ] . By sub- 
stituting the linear approximation for j?(y^) from equa- 
tion (3.2.10) into equation (3-2.1), a linearized approx- 
imation of the nonlinear differential equation results: 

Of' 


1 


Z (t) = 


Of 

3y 


/(t) +( F[y°(t)] - 


Since 


OF’ 

'Oy 


o 


and J° (t) are known functions of time, equa- 


'Sy 


y°(t) 


o 


( 3 . 2 . 12 ) 


o 


tion (3.2.12) is recognized as a vector linear differen- 
tial equation with time varying coefficients. The quan- 
tity in brackets is the driving function, a known function 


of time. 
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It is assumed that the solution of this linear- 
ized differential equation will give an approximate solu- 
tion to the nonlinear dif f erential equation subject to a 
given set of initia.l conditions. In particular,, suppose 
that the solution y^(t) found by integrating equations 
(3.2.12) is, in some sense, a better approximation to 
the solution of equation (3.2.1) than y°(t). If this is 
the case in general, an even better approximation for 
the solution of the nonlinear equation could be obtained 
by discarding y°(t) and integrating the linearized equa- 
tion again, this time regarding the new solution y (t) 
as the solution estimate. Evidently this process could 
be repeated any number of times. Iteration (n + 1) can 
be written as 


y 


1 n+1 


dF' 

3y 


y n+1 + 


h 


I(y n ) 


dP 

3y 


n 


y 


n 


(3-2.13) 


The superscripts here are iteration indices and do not 
represent exponentiation. 

The iteration represented in equation (3-2.13) 
is termed a generalized Newton- Raphs on iteration, since 
it may be obtained from a direct generalization of the 
Newton- Raphs on operator for finding roots of a scalar 
equation (McGill and Kenneth, 1964, p. 1762) . The 
question of convergence of this sequence has been pur- 
posely avoided, since even heuristic arguments are quite 
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complex. It is again noted that a sufficiency theorem 
for convergence does exist for systems which satisfy cer- 
tain stringent requirements on the set of differential 
equations and on y°(t) (McGill and Kenneth, 1963). 


3.2.2 Problems with Fixed Boundary Coordinates . 

Up until this point the boundary conditions of equation 
(3-2.3) have been ignored. The above arguments have, 
simply shown that a series of solutions to a linear 
initial value problem will, under appropriate circum- 
stances, converge to the solution of a nonlinear .initial 
value problem. But now, with linear differential equa- 
tions, the two-point boundary value problem is tractable. 
If the known final coordinates and the final time are 
fixed, the principle of superposition may be applied to 
equations (3-2.13) to satisfy the boundary conditions 
represented by equations (3-2.2) and (3-2.3)- For con- 
venience rewrite equation (3-2.13) for the general 
(n + 1) iteration as 


£ = Ay; + B 


where 


A = 


r^Fi 






(3-2.14) 


(3-2.15) 


(3-2.16) 
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At this point it is convenient to consider that the ele- 
ments of y (and., of course, F) have been ordered so that 
the p variables having known initial conditions are the 
first elements of y and the (2n - p) states having . 
unknown initial conditions are the last elements of y. 
This is in no way necessary. It Is simply notationally 
and computationally expedient. 

The homogeneous part of equation (3-2.14) is 

y = Ay (3.2.17) 

Consider integrating this equation forward (2n - p) times 

with the initial vectors chosen as follows: For each 

integration, each element of the initial vector is zero 

except for one. For the first integration, element 

th 

p + 1 is chosen to be 1. For the n integration, the 
(p + n)-th element is chosen to be 1. The reason for 
using these particular initial vectors to generate inde- 
pendent solutions to equation (3-2.17) will become 
apparent later. 

Let H be a matrix of resulting homogeneous solu- 
tions whose general element a^ is y^(t) generated with 

initial starting vector v defined above. In other 

k 

words, the k-th column of H is the solution vector y(t) 

corresponding to initial vector y . H has (2n - p) 

“°k 

columns and 2n rows. Since equation (3-2.14) is linear, 
the principle superposition applies, and a general 
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solution to (3.2.14) can be written as the sum of a par- 
ticular solution to the non- homogeneous equation and 
(2n - p) linearly independent solutions of the homogene- 
ous equation. In particular 

y(t) = H(t)ci + x(t) (3.2.18) 

where £ is a vector of arbitrary constants and where x(t:) 
is any particular solution of (3.2.14), which satisfies 
initial constraints of equation (3.2.2). In particular let 
x(t) = y (t), the most recent estimate for the optimal 
solution. At the final point this becomes 

y f = H(t f )o + x f (3.2.19) 

Equation (3-2.19) represents 2n equations in the (2n - p) 
arbitrary c's. However, p of these equations are of no 
interest since they involve variables which are not con- 
strained. Form the vectors y^ and x^ by deleting the p 
elements corresponding to states not involved in the 
final constraints. Form H' (t^) by deleting corresponding 
rows of Hp. H 1 is a (2n - p) by (2n - p) matrix. Then 
the pertinent equations may be written 

Zf = H’ (t f )e + Xf (3-2.20) 

The elements of c_ in equation (3-2.20) can now be selected 
to satisfy the required final boundary conditions. 
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_c = [H' (t f )] -1 [y^, - Xfl (3.2.21) 

The corresponding initial conditions are found by evalu- 
ating (3.2.18) at t and substituting jc from equation 
(3.2.21). Note that the first p rows of H(t ) are all 
zero, while the last (2n - p) rows and columns form an 
identity matrix. Thus, the vector _c does not affect 
known initial conditions and equation (3-2.18) becomes 
simply 

y io = c i + x io 1 = P+1. P+2, ... ,2n (3.2.22) 

where the c^ f s are given by equation (3.2.21). The 
advantage of choosing the initial vectors as described 
following equation (3.2.17) is now apparent. The y iQ 
values are referred to as the updated initial values, and 
the c^ values are referred to as the updates. 

If equation (3.2.14) is integrated forward, with 
the updated initial values from equation (3.2.22), the 
desired final endpoints will be approached. Under appro- 
priate circumstances, the resulting endpoints will be 
closer In a Euclidian sense to the desired final endpoint 
than that of the old solution. 

A complete iteration in solving the boundary value 
problem has now been made, and the entire process can be 
repeated again as follows: 
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(1) For a given set of initial values , repeatedly 
evaluate new values for y using equation. 
(3-2.13) until the solution of the linear dif- 
ferential equation (3-2.14) is a close approxi- 
mation to the solution of the nonlinear differ- 
ential equation (3-2.1). 

(2) Form the H matrix by integrating equation 
(3-2.17) with the appropriate initial vectors. 
Use the resulting final values to update the 
initial unknown values using equations (3-2.21) 
and (3-2.22) . 

(3) Repeat until the boundary values and the solu- 
tion are within specified tolerances. 

Note that step (1) may require 5 to 10 integra- 
tion iterations. Experience indicates that it is unneces 
sary to obtain an accurate approximation to the nonlinear 
differential equations before updating the initial values 
In practice , step (2) is performed after every integra- 
tion of equations (3-2.14) and (3-2.17)* except when two 
consecutive particular solutions differ considerably. 
Consequently, the convergence on the boundary conditions 
is simultaneous with the convergence on the solution. 

3-2.3 Variable Final Point Boundary . The method 
for variable endpoints is essentially the same as that 
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for fixed endpoints. Only a modification in obtaining 
the initial value updates c_ need be made. 

Two methods have been suggested for handling such 
problems with an unspecified final time. Long ( 1965 ) 
suggests a change in independent variable, t = as, where 
s is the new independent variable with a range of 
0 < s < 1. The variable a is appended to the system as 
a pseudo-state variable whose derivative with respect to 
s is zero. The method works well. However, a relatively 
complex term corresponding to the new state variable, a, 
must be added to each differential equation. A second 
method has been developed by Lewallen ( 1967 ) and is the 
method followed below. 

If the final time is free, there must be 2n + 1 
boundary conditions to completely specify the solution. 

If, as before, there are p initial coordinate constraints, 
there must be (2n + 1 - p) constraints involving the . 
final point. Let these be represented by 

t ± (yp = i = 1,2, ...,2n+l-p (3.2.23) 

or 

Jt(Xf ) = if 

where represents a constant vector of desired values 
for the final constraints. If the terminal constraints 
(3.2.23) are not satisfied, the difference between the 
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desired value j<* and the actual value jj_(x^) is called the 
terminal constraint dissatisfactions A#_: 

AJl = if - i(yp (3.2.24) 

The terminal constraint dissatisfactions can be linearly 
approximated by summing 

(1) the change in Jj_ due to small change in its argu- 
ments Ay^; assuming the final time fixed and 

(2) the change in j/_ due to a change in the final 
time At^ assuming the final state fixed. 


Symbolically this is 



' 8jj_ -j 


djj_ 



Ajj_ = | 


A Zf 

+ dt 

At 

If 

(3.2.25) 


An equation has already been obtained from which Ay^, can 
be expressed. Rewriting equation (3.2.20) gives 

= Xf ■“ 2Lf “ H ’( t f )£ (3-2.26) 

Substituting (3.2.26) into (3.2.25) yields 



' - 

d# 



A_Z - 

J 


At „ 
f 

(3.2.27) 


The values 
values of • 


for 

dy; 1 

dt , 


dj_ 

"dt 


can be computed directly from the 


using the differential equality 


di 

___ 


dy ! ' 

dt 

f 

__ d Zf 

dt 


f 


(3.2.28) 



d y ’ 

The values for the vector -rr— 


are available as. 


f 

a result of the last step taken in the numerical integra- 
tion of the particular solution. Substituting equation 
(3-2.28) into equation (3-2.27) and simplifying yields 


Ajj_ - 

In practice, 
matrix H 1 (t^) 


dy r 

~ 

H ( fc f + dt 

At 
f ■ 




the column vector 


dt 


(3-2.29) 
is adjoined to the 


and equation (3.2.29) is written 



' djf_ 

r 

dy" 

Ajr_ = 

_3yJ _ 

H’(t f ) 

dt 


(3-2.30) 


where 



The updates for the initial values, c_ and At^, can be 
solved for using equation (3-2.30) by a simple matrix 
multiplication and matrix inversion. 

3-3 -An Example of the Modified Newton-Raphson Method 

As an example of theory developed In the previous 
section, consider the solution of the two-point boundary 
value problem arising from the optimal orbital mechanics 
problem posed in Chapter 5- To apply the modified 
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Newton-Raphson methods the set of differential equations 
to be integrated must first be linearized. In terms of 
the theoretical development , the elements of matrix A and 
vector B, defined by equations (3.2.15) and (3*2.16), 
must be determined. For illustration, consider the u 
equation, equation (5-1-50) > which is repeated below for 
convenience . 


u 




(3.3.1) 


The elements of the vector y for this problem are 

m 

[u,g,R,\ ,1 ,\ j - According to equation (3-2.15) the 
partial derivatives of f^ with respect to the elements 
of y form the first row of the A matrix. These deriva- 
tives are listed below. 


An n — 


An 0 — 


A 


13 


A 


14 


3f 1 

_ 3f i 


3u 

Sf 1 

_ af l 

0y 2 

dg 


3f i 

^3 

" 6r 

3f l 

_ Sf l 

^4 

' 'K 


XX. 


M 


u 


3 


D 


JL 

R‘ 


2 cos § 


R 


3 


sin g 


1 ig_ p-3/2 

M 2 v 
u 


(3-3-2) 

(3-3-3) 

(3-3-4) 


(3-3-5) 
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B vector are exhibited in FORTRAN in subroutine DIFFL of 


Appendix B. 

Next, the boundary conditions must be linearized, 

dt 

To do this the elements of the matrix 




must be 


determined. For the optimal orbital mechanics problem. 


[u f ,g f , R f ,X uf ,\ gf A rf ] 


T 


since all of the elements 


of y^ are involved in the constraints. Following the 
notation used in Appendix B, this matrix is called HPAR. 
The constraint equation (5-1*29) is repeated below as an 
example . 


u 


f 


1_ 

R, 


C 


E 


(3-3-11) 


The elements of the first row of HPAR are given below. 

d\J/ di!/-, 

HPAR 11 = -^~r = = u f 

diK di|/ 

HPAR no = - 0 

12 dy 2f 

df -i dilf -i 

HPAR 13 " ByJ = ^ = ^2 

HPAR^ = HPAR 15 = HPAR l6 = 0 

All the elements of the HPAR matrix are exhibited in 
FORTRAN in subroutine BOUND of Appendix B. Once the pre- 
liminary work of linearization has been performed, the 
generalized subroutines exhibited in Appendix B can be 
used directly. 


(3-3-12) 

(3-3-13) 

(3-3-14) 

(3-3-15) 



Subroutine QUASI controls the generation of a 
solution to a two-point boundary value problem. Although 
the listing of this subroutine in Appendix B is intended 
to be self-explanatory, a general description of the com- 
putational procedure is in order. Each generalized 
Newton- Raphs on iteration requires the generation of a 
particular solution and (2n - p) homogeneous solutions to 
equation (3-2.14). Since each solution has 2n independent 
variables , one iteration requires the integration of 
2n(2n - p + 1) first order differential equations. Each 
of these solutions could be generated separately. How- 
ever, if rapid access memory is limited, this method is 
computationally inefficient since matrix A(t) must be 
evaluated (2n - p + 1) times. Instead, all (2n - p + 1) 
solutions are generated simultaneously. At any given 
time in the integration, matrix A need only be evaluated 
once in order to integrate all solutions forward one 
step. This method also saves considerable time in com- 
putational overhead. 

Before solving the two-point boundary value prob- 
lem, an estimate of the initial solution is stored in 
matrix XOLD at discrete times differing by a constant 
increment A. XOLD is defined so that XOLD(i,l) = 
y ± (t 0 ) and XOLD(i,j) = y ± (t 0 + (j-l)A). 

A summary of the operation of subroutine QUASI is 


given below. 
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(1) Set the elements of matrix A and matrix B which 
are not functions of time to their initial val- 
ues (call DIFFl). 

(2) Set up initial vector X of length 2n(2n - p + 1) 
to he integrated. The elements of X correspond 
to the initial vectors for the particular solu- 
tion and (2n - p) homogeneous solutions needed. 

(3) Integrate all equations to the final time esti- 
mate. After each integration step, store the 
new particular solution in XOLD. The old par- 
ticular solution, which has just been used to 
evaluate A and B, is of no further use and is 
destroyed . 

(4) Evaluate the updates to the unknown initial 
values and to the final time estimate (call 
BOUND) . 

(5) Change the unknown initial conditions and final 
time by a fraction of the nominal updates . The 
fraction used depends upon certain convergence 
and stability criterion. 

(6) Determine if the solution has converged to 
within desired tolerances. If not, return to 
step 3- 

The actual integration in step 3 is performed by 
subroutine RUNKUT which calls subroutine DLSUB to 
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evaluate the right hand side of the first order differ- 
ential equations being integrated. DLSUB in turn calls 
subroutine DIFFL, which evaluates matrices A and B using 
the old particular solution stored in XOLD. Using A and 
B, subroutine DLSUB then evaluates the right hand sides 
of the particular and all homogeneous differential equa- 
tions and returns these results to RUNKUT. Subroutine 
RUNKUT integrates forward one step and returns control 
subroutine QUASI, which stores the new solution in XOLD 
and determines if the final time has been encountered. 
This procedure is repeated until the final time is 
encountered. The details of these programs are described 
in the FORTRAN listings of Appendix B. 



CHAPTER 4 


THRUSTING HARMONIC OSCILLATOR 

This chapter deals with the occurrence of mul- 
tiple stationary solutions for problems having bounded 
control and periodic solutions. Such multiple solutions 
have been observed in a nonlinear orbital mechanics 
problem discussed in Chapter 5- In order to gain some 
insight into the latter problem, it is desirable to 
investigate first a simple linear problem which exhibits 
similar multiple solutions. 



Fig. 4.1 Thrusting Harmonic Oscillator 
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Consider a mass m which is connected vertically 
to an inertial reference by a spring and dash pot as 
shown in Figure 4.1. The spring constant per unit mass 
is represented by k, and the coefficient of damping per 
unit mass is represented by c. The mass is capable of 
producing a bounded thrust per unit mass of T in the 
upward direction only. The mass is assumed constant. 
The equations of motion for this thrusting harmonic 
oscillator are 


where 


v = T — cv - kx 
x = v 


0 < T < T 


max 


k > 0 , c > 0 


(4.1.1) 

(4.1.2) 

(4.1.3) 


x is the position of the mass measured in an 
upward direction from its equilibrium posi- 
tion 

v is the mass velocity. 


The problem is to transfer the mass from rest 
to a given final height h while minimizing the control 
effort, 

h 

J = / Tat (4. 1.4) 

t 


O 
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The boundary conditions (algebraic equations of con- 
straint) for this problem are 

x^ = 0 v = 0 t = 0 (4.1.5) 

o o o v ' 

= h v^ is free t^ is free (4.1.6) 

The Hamiltonian and G functions may be written as 
H = X (P - cv - kx) -f \ v - F 

v .X 

= (X y - 1)F - X y (cv + kx) + X x v (4.1.7) 

G = h 1 (x f - h) ' (4.1.8) 

The initial point constraints have not been included in 
the sum G above, since doing so provides no useful infor- 
mation. Since H is not an explicit function of time, 
the first integral implies that H is equal to a constant. 
The adjoint variable Euler- Lagrange equations are 

X = X c - X (4.1.9) 

and 

X x = X y k (4.1.10) 

The control switching function is seen from equation 
(4.1.7) to be 


1 
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Therefore , the Maximum Principle (see section 1.3-5) 
yields the following control algorithm: 


T = 0 


0 < T < T 


max 


T = T 


max 


X v " 1 < 0 
X v ” 1 3 ° 

X V - 1 > 0 


(4.1.11) 


For singular control, that is, intermediate thrust, 

X must equal I, and all the time derivatives of X must 
v ^ v 

be 0. If this is the case, equation (4.1.9) requires 
X = c; therefore, X is zero and equation (4.1.10) 
requires k = 0. This contradicts the assumptions of the 
problem and singular control is ruled out. The control 
must be bang-bang. 

Solving equations (4.1.9) and (4.1.10) for X y 


gives 

—t 

X y = e 2 (A sin bt + B cos bt), (4.1.12) 

if 

c 2 - 4k < 0 


where 

b = ~(4k-c 2 ) 1//2 


(4.1.13) 


2 

The case of c - 4k _> 0 is not of interest since the 
physical system does not exhibit periodic oscillations 


under these circumstances. 
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The possibility of null thrust 
T(0 + ) = 0 


(4.1.14) 


for the initial control will not be considered since the 
state would remain unchanged. Starting with maximum 
thrust it is seen from equation (4.1.7) with H(0) =0 
(to be shown) that 




(4.1.15) 


In order to maintain T = T „ at the next instant of 

max 

time j 


X (0) = \ vo > a, positive constant (4.1.16) 

The function X y is shown in Figure 4.2 for a 
typical set of parameters. Note that the envelope of max- 
imum points is asymptotically increasing. The correspond- 
ing control is shown in Figure 4.2b. Notice that T = T max 

for X >1, and T = 0 for X <1. The control is seen to 
v v 

be a series of thrusting intervals separated by coasting 
or null thrust intervals. Each thrusting interval is 
longer than the previous one. However, since the period 
P of X y is fixed, no thrusting interval can be longer 
than the P/2. The transversality conditions require that 


X 


vf 


0 


(4.1.17) 


H f = 0. 


(4.1.18) 
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The former condition is the cutoff condition and the lat- 
ter “condition requires that H = 0 throughout. If these 
conditions are satisfied , equation (4.1.7) gives 

X xf v f - F f = 0 (4.1.19) 

Since X v ^ = 0 by equation (4.1.17) > it can be seen from 
equation (4.1.11) that 

F f = 0 (4.1.20) 

Since Is always negative (see Figure 4.2), equations 

(4.1.17) and (4.1.9) determine that X xf ^ 0. Thus, 
equations (4.1.19) and (4.1.20) require that 

v f = 0 (4.1.21) 

Although analytic solutions can be obtained, the 
evaluation of the boundary conditions becomes tedious for 
solutions with multiple thrusting intervals. The nature 
of the solutions Is easily observed with the aid of 
numerical integration. A parametric flooding technique 
(see section 3-1*1) is practical, since solution curves 
are members of a one parameter family of curves with 
parameter X XQ , or equivalently, x vo - That is, all pos- 
sible solutions can be generated by continuously varying 

X from 0 to + . 

vo 

With a value chosen for X yo ^ differential equa- 
tions (4.1.1), (4.1.2), (4.1.9) and (4.1.10) are 



Switching 
Value 1 


Cutoff 0 
Value 


T 

max 


0 



t 


t 


3 Thrust 

Pinal 

Point 


t 


Fig. 4.2 The Cutoff Function., Thrust Magnitude 
and Height 
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integrated until one of the zeros of the cutoff function, 
equation (4.1.17);- is encountered. As described in sec- 
tion 1.3.4, the zeros of the cutoff function represent 
points satisfying the transversality conditions . 

Since v is obviously a periodic function for 
c - 4k < 0, equation (4.1.21) implies that the cutoff 
function is satisfied during each period of oscillation. 
Since a thrusting interval occurs during each period, 
solutions exist with any number of thrusting periods. 

The existence of multiple stationary solutions 
is easy to observe. Limits exist on the range of posi- 
tive heights that may be reached by a solution with any 
specified number of thrusting intervals. For a solution 
with two thrusting intervals, the minimum specified 
height which can be optimally reached may be found by 
setting A_ vo = € and numerically integrating until the 
third zero of the cutoff function \ occurs. It can be 
seen from Figure 4.2 that as e approaches 0, the length of 
the first thrusting interval approaches 0. On the other 
hand, as A. vq approaches + °>, the length of the initial 
thrusting period approaches P/2, the maximum allowable. 

In this way the range of achievable positive heights can 
be computed. 

Figure 4.3 shows the range of heights which can 
be reached for solutions with up to seven thrusting 
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intervals. Note that al3. heights between .13 ft. and 
some maximum attainable height can be reached in more than 
one way. Heights between .41 and 1-99 ft. can be reached 
by solutions with one, two, or three thrusting periods. 
Each of these solutions is a stationary solution satisfy- 
ing the first three necessary conditions of section 1.2. 
Here is a clear case of multiple stationary solutions . 

The number of stationary solutions to a given height 
increases as the given height increases. From Figure 
4.3 it can be seen that a height of 5-77 ft- can be 
reached with from three to seven thrusting intervals. 

Although each of these solutions satisfies neces- 
sary conditions for optimality. Figure 4.4 shows that 
the effort required to achieve a given height is differ- 
ent for each mode of transfer. Only one global optimum 
exists in each case. To reach a height of 5-77 ft. 
requires an effort of 9 -42 lb.p-sec. using three thrusting 
intervals. By using the optimal four thrusting intervals, 
the effort can be reduced to 8.36 lb f -sec. Five thrust- 
ing intervals requires an effort of 8.80 lb^-sec. 

Although not verified analytically, it is evident 
that each solution to a given height satisfies the end- 
point sufficiency condition developed in Chapter 2. 

That is, each solution represents a local minimum of 
effort with respect to small variations in the final ■ 
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velocity. This conclusion may be drawn from energy con- 
siderations. From an energy standpoints the oscillator 
will gain the maximum amount of energy by centering each 
thrusting period about the maximum positive oscillator 
velocity. Therefore, if the length of each thrusting 
period is fixed at the optimum value, but the time of 
thrust initiation is changed slightly, the oscillator 
will gain less than optimum energy from each thrust and 
consequently attain a less than optimum height. 

Without further sufficiency conditions, no cri- 
teria exist for selecting the global minimum for this prob- 
lem. There seems to be little hope for the development 
of global sufficiency conditions from the calculus of 
variations in its current formulation. The reason for 
this is apparent. The necessary and sufficient conditions 
that have been developed to date are effective only in 
determining trajectories which are local optimums. That 
is, only small variations in the state variable trajec- 
tory are considered. It is clear that the state trajec- 
tory for a two thrust control program is not a small 
variation from the state trajectory for a single thrust 
program. 

Referring to similar multiple thrusting impulses 
involved in optimal orbital transfers, D. F. Lawden (1963, 
pp. 112-111) concludes the following: 
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. . . the conditions we have found for an 
optimal trajectory serve only to identify 
maneuvers which are optimal relative to small 
variations of the thrust program. This implies 
that a number of such optimal trajectories may 
be available in any particular case and we have 
no criterion., other than direct comparison of 
the characteristic velocities [performance index] 
of each; for deciding which represents the ab- 
solute optimal mode of transfer .... 

. . . modes employing two or more impulses and 
satisfying all of the [necessary] conditions may 
also be found and these , also, will represent 
relative optima. It is clear that the theory, in 
its present state, requires that we should first 
-choose the order in which the various thrust 
phases shall succeed one another, and our [neces- 
sary] conditions will then select the relative 
optima from the class of all programs for which 
the thrust sequence follows the chosen pattern. 

No criteria is known for which the optimal pat- 
tern can be ascertained beforehand. 



CHAPTER 5 


AN OPTIMAL ORBITAL MECHANICS PROBLEM 

This chapter considers the problem of transfer- 
ring a low thrust space vehicle from an initial circular 
orbit to any coplanar elliptic orbit of given energy and 
angular momentum. The problem is to determine the thrust 
program so as to minimize 

t f 

J = j* T(t) dt 
t 

o 

where T(t) is the magnitude of the thrust. A typical 
transfer is shown in Figure 5 * 1 * The final argument of 
periapsis the final true anomoly cp^? the final range 

angle f^, and the final time t^ are all unspecified and 
considered free. For such transfers from an initial 
circular orbit , the orientation of the final orbit is of 
no concern, since once is known, any desired orienta- 
tion can be optimally achieved by simply coasting in the 
initial circular orbit for an appropriate period of time. 

Many investigators have turned their attention 
to optimal orbital transfers in the past two decades . 

Bell (1968) gives an excellent survey of problems which 
have been investigated, including a bibliography with 
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160 entries. Recent investigations concerning the mini- 
mum fuel problem have followed two distinct approaches. 

In the first approach, simplifying assumptions are made 
which allow analytical treatment of some problems . 

Lawden ( 1963 ) , Robbins (1965), and others have 
analytically treated minimum fuel transfers using an 
impulsive approximation. McIntyre and Crocco (1966; 

1967) have examined minimum fuel transfers between circu- 
lar orbits under the condition that the separation of 
the orbits is small and have found analytic approximations 
exhibiting multiple thrusting periods. 

Problems concerning minimum fuel orbital trans- 
fers with low thrust have usually been investigated 
using a computational approach. Numerical studies of 
this problem have been made by McCue (I967), Handelsman 
(1966), and McCue and Bender (1965) . Both analytic and 
numeric investigations have contributed substantially to 
the understanding of the problem of minimum fuel orbital 
transfers . 

In the first section of this chapter, this partic- 
ular problem is set up and the optimization conditions 
derived. An instance of multiple stationary solutions 
not noted by previous investigators is described in sec- 
tion 5 - 2 . In section 5-3 the endpoint sufficiency condi- 
tion for this problem Is described and the computational 
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procedure for implementing it is discussed. In section 
5-4; several cases of multiple solutions are illustrated 
with trajectories in semi-major-axis-eccentricity space. 
The final section describes the general nature of the 
low thrust orbital transfers which were investigated. 


5.1 Analysis 

Using a natural coordinate system, the differen- 
tial equations governing the dynamics of a space vehicle 
shown in Figure 5-1 are 


v 

,.g 

r 

m 

f 


* 

T GM 

- cos x g— sm g 


T . GM v 

— sm x - — 7T- cos g + — cos g 
mv 2 r & 

r v 


v sm g 
T 

v _ 


v 


cos g 


(5.1.1) 

(5.1.2) 

(5.1.3) 

(5.1.4) 
(5-1.5) 


where 

v = velocity magnitude 

g = flight path angle (the angle from the local 
horizontal to the velocity vector) 
r - radial distance from the gravitational force 
center to the vehicle 
m = vehicle mass 
f = range angle (f(t ) =0) 
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x - thrust vector control angle (the angle from 

the velocity vector to the thrust vector) 

T = thrust magnitude 

G = universal gravitational constant 
* 

M = mass of the gravitational force center ; 
constant 

v g = engine exhaust velocity; constant 

In deriving the above equations , only planar 
motion about a homogeneous , spherically symmetric central 
body has been considered. The vehicle has been assumed 
to be a point mass with mass much less than that of the 
gravitational force center. In addition, it is assumed 
that the thrust magnitude is bounded by 


0 < T < T 


max 


( 5 . 1 . 6 ) 


It is useful to nondimens ionalize equations 
(5.1.1) - (5*1.6) in order to compare orbital transfers 
about one central body with those about another. The 
non-dimensionalization is based upon the energy of the 
initial orbit E q and variables associated with a circu- 
lar orbit having the same energy. 

For this purpose define the following constants 
GM*m 

■p — 

c 2E 

o 


(5.1.7) 
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(5.1.8) 



(5-1.9) 

( 5 . 1 . 10 ) 

(5.1.11) 


h = R v m (5.1.12) 

c c c o ' ' 

where 

m = initial vehicle mass 
o 

a = acceleration constant 
c 

h c = angular momentum constant 

With these constants, the following nondimensional vari- 
ables are defined 

R = — — = nondimensional radius (5-1.13) 

c 

u = = nondimensional velocity (5.1.14) 

v c 



nondimensional time 


(5.1.15) 


M = = nondimensional mass ( 5 .I.I 6 ). 

o 


T 

P = — - — = nondimensional thrust 
m a 
o c 


(5-1.17) 



E = r=r- = nondimens ional energy 

iii 

c 


h = = nondimens ional angular momentum (5- 1.19) 

c 


Substituting equations (5.1*13) - (5-1*19) into the sys- 
tem differential equations yields the following nondimen- 
sional system equations: 


jp 2 . 

u' = cos x - -75- sin g 


( 5 - 1 . 20 ) 


P 1 

g‘ = sin x - — g— cos g + ^ cos g 


( 5 . 1 . 21 ) 


R 1 = u sin g 


( 5 . 1 . 22 ) 


(5-1.23) 


f ' = ^ cos g 


(5-1-24) 


where the prime indicates differentiation with respect to 
the nondimens ional time t • 

The problem to be investigated then is that of 


minimizing 


J = J F dr 


(5-1-25) 


subject to differential constraints (5.1.20) - (5.1.24) 
and endpoint constraints 



1 = 0 


115 

( 5 - 1 . 26 ) 



(5-1.27) 

(5-1.28) 

(5-1.29) 

(5-1.30) 


where C,-, and C, are constants . 

E h 

Constraints (5-1.26) - (5.1.28) specify an initial circu- 
lar orbit. Equation (5-1.29) specifies that the energy 
of the final orbit is to be C R and equation (5.1-30) 
specifies that the angular momentum of the final orbit is 
to be C R . The final energy and angular momentum specify 
the shape of final orbits but not the orientation. 

To determine the necessary conditions for a mini- 
mum fuel orbital transfer, the H and G functions, equa- 
tions (1.2.7) and (1.2.8), are first formed: 


F 1 

H = l u [m COS X 2 sin 


, , ,P . 1 

+ X g ^ Mu x - cos 

& mu R ^ u 


+ COS g] (5-1-31) 




1 16 


Equation (5.1.24) need not be included in H, since f is 
an auxiliary variable not directly coupled to the other 
equations . The G function is 

G = |i x ~ C E ] + P 2 t R f u f cos Sf ~ C h^ (5-1.32) 

Constraints (5-1.26) - (5-1.28) have not been included in 
the G function, since doing so would yield no useful infor- 
mation. Forming the adjoint variable Euler- Lagrange 
differential equations,, equations (1.2.10), yield 

F 1 1 

X' = [ — 9 Sin X 9-0 cos g - p COS g] . (5.1.33) 

u s Mu R u K 


- X r sin g 


xi = 


c os g 
u R 2 


r_4 

1 .2 


u-, 

- sm g 


R u 


(5-1.34) 


X u cos 
r 


K - 


2\ 


u 


R 


3 


sin g - X ^p] cos g 

g R J u R 


(5.1.35) 


X' = cos x + L — sin X (5-1.36) 

m U g Mu 


The control variable Euler- Lagrange equation, 
equation (1.2.11) for control x is 


§ [ - X u sin x + ^ cos x] - 0 (5-1.37) 
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For F £ 0 equation (5-1-37) implies 
sin x = 


and 




(5-1-38) 


(5-1-39) 


Since the second control variable F is bounded 
and appears linearly in H, "bang-bang" control may be 
expected. The possibility of singular control described 
by Lawden (1963, pp. 86-94) has been determined to be 
nonoptimal, in general, by Kelly, Kopp, and Moyer (1967, 
pp. 92-100). Therefore, the case of the switching func- 
tion S = 0 is not considered here. By substituting 
equations (5-1-38) and (5.1.39) into equation (5-1-31) 
and simplifying the term multiplying F, the thrust 
switching function may be written 


S 


1 1 

,2 p 

f x 

g' 

+ X - M 

' 1 

l u i 

I u 

u 

1 e 


(5-1-40) 


The thrust control may be written as 


F = F„ 


max 


S > 0 


F - 0 


S < 0 


(5-1-41) 
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The transversal! ty equations , 
15), yield 

equations (1.2.12) 

^l u f + ^2 R f cos g f + x uf = 0 

(5-1-42) 

|a 2 R f u f sin g f + X f = 0 

(5-1-43) 

+ q 2 u f cos g f + X rf = 0 
R f 

(5-1-44) 

Vf = 0 

(5-1-45) 

0 

1! 

K 1-1 

(5.1.46) 


Eliminating and p.g from equations (5.1.42) - (5-1.44) 
yields 


sin g f = 0 (5.1.47) 

(5-1-31 ) , it is observed 
that condition (5-1-46) becomes 


'gf 


u. 


1 


R f u f R f 2 / 


COS 


+ (l f»U_p - p) 


Rr 


Using equations (5-1-47) and 



- *r (1 + 


u 


e 


= 0 


(5-1-48) 
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Equations ( 5 . 1.45) j (5- 1.47) and ( 5 . 1.48) constitute a 
set of necessary conditions which the endpoints must 
satisfy. Comparing equation (5.1.40) to equation 
. (5-1.48) , it is evident that points which qualify as 
endpoints also qualify as switching points for the thrust 


magnitude control. 

Imposing condition (5.1.45) on equation (5-1.48) 
yields a cutoff function. 



(5.1.49) 


The zeros of this function determine points at which 
numerical integration may be terminated. Substituting 
the optimal thrust direction control from equations 
( 5 . 1 . 38 ) and (5.1.39) into the nondimensionalized state 
variable differential equations > (5.1.20) - ( 5 - 1 . 23 ) and 
adjoint variable differential equations., (5- 1-33) - 
( 5 .I. 36 ) gives 

FX -1 

u' = — =rjo ? sln S (5.1.50) 

MD 7 R 




cos g 


(5.1.51) 


R 1 = u sin g 

F 
u 

e 


(5.1.52) 


M 1 


(5-1.53) 
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K - \ 




i__ + i 

2 2 + R 


R u 


cos g 


X sin g 


X cos g 

i _ u 


R' 


2 


1 

R 2 u 


u 

R 


sin g - X r u cos g 


2X u 

1 ' = -q- sin g - X 

I R-J i 


2 u 


R 3 u 


R 


cos g 


( 5 - 1 . 54 ) 

(5.1.55) 

(5.1.56) 


X 


I 

m 




(5.1.57) 


where 


The integration of this set of equations subject 
to equations (5. 1.26) - (5.I.30) and (5-1.42) - (5.1.46) 
yields the desired optimal trajectories. 

5-2 Existence of Multiple Stationary Solutions 

The problem posed in the previous section per- 
mits the existence of multiple thrust solutions . Mason 
(1967, pp. 49-70) has shown the existence of two solu- 
tions to a problem as posed in section 5-1* one solu- 
tion is characterized by a single continuous thrust to 
the desired endpoint; a second solution, composed of two 
thrusting periods separated by a null-thrust coasting 
arc, reaches the same endpoint with a different 
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performance index. Both solutions satisfy the necessary 
conditions of the calculus of variations . This type of 
multiple stationary solution is similar to the type of 
multiple stationary solutions exhibited by the thrusting 
harmonic oscillator in Chapter 4. Because of this simi- 
larity,, it is probable that both the single thrust and 
the thrust-coast-thrust transfers satisfy the endpoint 
sufficiency condition described in Chapter 2. Verifica- 
tion of this hypothesis is left to future investigations. 
Several other investigators have also observed optimal 
orbital transfers with multiple thrusting periods 
(McCue, 1967i Handelsman. 1966; McIntyre and Croce o, 
■ 1967 ) . 

The investigation pursued in this chapter will 
instead be limited to orbital transfers with a single 
thrusting period. That is, thrust magnitude P is not 
considered to be a control, but will be of fixed magni- 
tude for the entire duration of the transfer, i.e., 


< t < t, 


(5.2.1) 


If P is not a control, the M equation, equation (5.1.53) 
is no longer coupled to the system of differential equa- 
tions represented by equations (5.I.50) - (5.1. 57). 
Instead, M is a known function of time independent of 
u, g, R, and x. 
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M = (t - t Q ) (5.2.2) 

e 

Therefore the M equation need not he adjoined to the sys- 
tem through the Lagrange multiplier All equations 

in the previous section are explicitly independent of 
time and will remain valid if is set to zero and M is 
evaluated using equation (5-2.2). 

Optimal trajectories satisfying equation (5.2.1) 
are determined hy the six differential equations (5.1.50)- 
(5.1.52), and (5.1.54) - (5 .I. 56) subject to the seven 
boundary conditions represented by equations (5. 1.26) - 
(5.1.30), (5.1.47) and (5-1.48). 

Throughout this chapter, orbits will be displayed 
as points in a plane having coordinates of nondimensional 
semimajor axis A and eccentricity e. Except for the 
sense of the angular momentum vector, there is a unique 
mapping between the variables E and h, whose final val- 
ues have been specified and the variables A and e. The 
relationships are given below 

A = - |g (5.2.6) 

e 2 = 1 - 2 Eh 2 (5-2.7) 

Here A is the nondimensional semimajor axis. The non- 
dimensional variables E and h are defined by equations 
(5.1.18) and (5.1.19). 
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Before a solution to the boundary value problem 
was attempted^ the method of flooding described in sec- 
tion 3-1 was used to generate a locus of optimal end- 
points in A-e space. To give physical significance to 
the required initial guesses, the following equations 
were derived for a fixed mass vehicle which relate the 
initial unknown Lagrange multipliers to the initial value 
of the control angle x q and its derivative : 


uo 


L go 


cos x 


o 


u sin x 
o o 


1 


TO 


u 


"cos (g^rr;) 


o • 


F_ 


' max 


sin x - u x 
o o 


(5-2.3) 

(5-2.4) 

(5-2.5) 


u cos g o 

+ ■£- sin x Q sin (g Q + x Q ) 1 


R. 


o 


The first two of these relations result from the control 
Euler -Lagrange equation (5- 1-37) , the initial conditions 
and the fact that H q = 0. The last equation results 
from the implicit differentiation of the control variable 
Euler-Lagrange equation and subsequent elimination of the 
resulting time derivatives using the state and adjoint 
differential equations. 

Figure 5-2 displays the locus of optimal end- 


points reached by initiating integration with x Q = 4 


o 
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and varying values of x^. The coordinates used in the 
figure are semimajor axis and eccentricity. Each point 
on the locus was determined by terminating integration 
when the first zero of the cutoff function (5.1,49) was 
encountered. As x^ was increased from a. given initial 
value , the final point of the trajectory moved continu- 
ously forming the locus curve OABC. With a further 
increase in x ^ the locus was generated from C to D pass- 
ing a second time through point A. A still further 
increase x^ resulted in a jump from D to E and then a 
continuous return to 0 along curve EO. The fact 'that 
the locus of optimal endpoints intersects itself at A 
indicates that two solutions exist to point A. 

The multiplicity of solutions is compounded if a 
second locus of optimal endpoint for x Q = 16 ° is super- 
imposed on the one for x q = 4 ° as shown in Figure 5 • 3 • 
The movement of the optimal final point along the x Q = 
16 ° locus for increasing x^ is similar to that described 
for the x q = 4 ° case. The l6° locus not only intersects 
itself , as did the 4 ° locus , but also intersects the 4 ° 
locus in three places. Since similar loci exist for all 
intermediate angles , 4 ° _< x q _< 16 °, it is evident that 
a complete region of endpoint space exists which can be 
reached by two distinct "optimal" trajectories. 
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5-3 Computational Procedures 

The endpoint sufficiency test developed in Chap- 
ter 2 was applied to the case of multiple stationary 
solutions described in the previous section. Since the 
nonlinear differential equations describing the orbital 
transfer problem can not be integrated analytically, 
resort must be made to the numerical algorithm prescribed 
for such cases by section 2 . 6 . This numerical algorithm 
assumes the user has the capability of solving two-point 
boundary value problems. Chapter 3 has set forth a 
method for solving such problems, complete with detailed 
references to allied computer programs listed in Appen- 
dix B. It remains for this section to describe in detail 
programs designed to implement the numerical algorithm 
proposed in section 2.6. In addition, the overall 
computational scheme for determining minimum fuel orbi- 
tal transfers will be outlined, thus completing the 
description of programs contained in Appendix B. 

5-3-1 The Endpoint Sufficiency Condition . The 
endpoint sufficiency condition is given by equation 
(2.3-5)- In this case the general endpoint is repre- 
sented by the eight variables ( u 0 , g , R 0 ,T 0 ^u f ,g f , R f ,r f ) - 
The endpoint constraints are represented by equations 
( 5 - 1 - 26 ) - ( 5 . 1 . 30 ), and the following constraint on 


the initial time : 
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(5.3.1) 


T 


O 


= 0 


Since there are six equations involving the eight end- 
point variables., there are two degrees of freedom. That 
is, two endpoint variables involved in constraints which 
relate two or more endpoint variables (nontrivial con- 
straints) may be independently specified, and the other 
six variables will then be given by the equations of 
constraint. Constraints which fix an endpoint coordinate 
simply result in a row of zeros in the Q matrix of equa- 
tion (2.3.5). Therefore, only the two nontrivial con- 
straints, equations (5. 1.29) and (5.I.30), and the non- 
fixed' coordinates u f , g^, r^ and will be considered 
in the following analysis. 

Any two of these variables can be selected as 
the independent variables. For convenience, r f and 
are chosen as the independent variables . 

In terms of the notation of section 2.3., the 
pertinent relations are 


V = 

Rf 

w = 

U f 


*f 


g f 


1 2 
2 u f 

R^u 


1 _ 

Re 


'E 


: f u f cos g - C 


h 


g f 

R f 

T f 


\ 

>(5.3.2) 


/ 
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However , as was the case in the example of section 2.5; 

J* need not be considered to be a function of since 
actual integration of the Euler- Lagrange equations pro- 
vides a relationship between x^, and the other final point 
variables. Hence v and r can be written more simply as 


u 


v - R , 


r = 


(5-3.3) 


R, 


Following the procedure described in section 2.3 yields 

1 

2 ! 

(5-3.4) 


cH 

cFv 


R, 


u f cos g f 


cH 

cFw 




R f cos g f 


0 


R f u f sin g f 


(5.3.5) 


and finally 


n 


u f R f 


cos g f .(l - u f ^R f ) 


u 


f R f sin g f 


(5.3.6) 


1 
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Following the form of equation (2.3-10) the matrix 


d 2 J*~! 

■ v- ' -v can 
o_roi^J 

matrices 


be written as the sum of the following two 


~ ~\2 
o J* 


d 2 G 

-4- l 

' d A^~ 

drdr 


drdr 

j 

dr 


(5-3.7) 


where A f = (X uf , X gf , X rf ) 

The above form will always result for problems 
with fixed initial coordinates which can be written so 
that J* is not a function of t f . The first term in equa- 
tion (5-3-7) can be determined from equation ( 5 . 1 . 32 ) 
yielding. 


a 2 c, 

d r 6 r 


^1 

P 2 R f sin g f 
P 2 cos g f 


q 2 R f sin g f 
ia 2 R f u f cos g^ 
p 2 u sin g f 


P 2 cos g f 

PgU sin g f 

2 l l l 

R 3 

(5-3-8) 


The variables p-^ and p 2 in the above expression are deter- 
mined from equations (5.1.42) and (5-1-43) giving 


»1 


p 2 [ X gf 
f [ R f tan g f 


+ X 


rf 


(5-3.9) 
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A gf 

R f u f sin g f 


(5-3.10) 


The second term of equation (5-3.7) is determined by com- 
putational procedure explained in connection with equa- 
tion ( 2 . 6 . 9 ). In this case represents one of the 
elements of A_ . 

The computational procedure is as follows: 

(1) The endpoint sufficiency test is initiated by 
a call for subroutine FOCAL. At this time 
matrix XOLD contains the time history of. a 
trajectory to the endpoint which is to be 
tested. 

(2) Subroutine FOCAL immediately calls subroutine 

CONSTRT. This subroutine used the nominal 

final point values r* from XOLD to evaluate 

— f r ^ 2 i 

. d G 

Q from equation (5-3 -o) and fnom equa- 

tions (5.3.8), (5.3.9) and ( 5 . 3 .IO). These 
values are returned to subroutine FOCAL. 


(3) 

( 4 ) 


Subroutine FOCAL then evaluates the matrix 
d — f in accordance with equation (2.6.9). 


Subroutine FOCAL then computes 


a 2 j* 


equation ( 5 - 3 - 7 ) and determines if it is 


positive definite. If the matrix is positive 
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definite , the endpoint being tested represents 
a local minimum. 


To evaluate 


in step 3 requires the solution of 


three two-point boundary value problems. For example, 

5 A 

to determine , a new final point is determined by 

making a small change A in u|. An optimal path to this 
new endpoint having fixed final state coordinates must 
now be determined. The original boundary value problem 
involved the variable boundaries of equations (5.1.29) 
and (5.1.30) and the resulting transversality equations 
(5-1. 47) and (5. 1.-48). For the new problem the final 
point constraints are 


u f = Cl 


(5.3.12) 


R f = c 3 (5-3.13) 

With these constraints the only useful relation given by 
the transversality conditions is 

H f = 0 (5-3.14) 

Both the original two-point boundary value prob- 
lem and the new problem can be solved using the same 
subroutine QUASI. This is due to the fact that QUASI 
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does not directly determine the initial value updates 
from the boundary conditions. Instead, QUASI calls a 
separate subroutine to handle this job. Since this sep- 
arate subroutine is an argument of subroutine QUASI, 
subroutine QUASI can be executed using any one of several 
different subroutines to handle differing sets of bound- 
ary conditions . The subroutine which handles the orig- 
inal variable final point problem is named BOUND. The 
subroutine which handles the fixed boundary values of 
equations (5. 3- 11) - (5 •3*1^-) is named FIXED. 

In the example, an optimal path to the new end- 
point is determined by calling QUASI with subroutine 
FIXED as an argument and with c-^ = u^ + A, c^ = g* and 
c^ = R|;. QUASI returns the new solution to FOCAL as a 
time history contained in matrix XOLD. Using the final 
values of the Lagrange multipliers from XOLD, the par- 

fl 

tial derivatives can be computed directly using 

L° u f J rSA ■> 

equation (2.6.9). The other two rows of ^ — J are com- 
puted in similar fashion. Subroutines FOCAL, CONSTRT, 
FIXED and BOUND are listed in Appendix B. 



solutions to specific two-point boundary value problems 
exhibited in the next two sections were found following 
the general computational outline given below. 
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(1) Generate an initial solution estimate by inte- 
grating forward the nonlinear equations using 
the known initial values and estimates for 
unknown initial values. The purpose of this 
integration is to determine an estimate of the 
final time. 

(2) With the final time known, it is possible to 
regenerate the Initial solution by integration 
with fixed step length. The solution is simul- 
taneously stored in the matrix XOLD at fixed 
equal increments in time, as required by sub- 
routine QUASI. 

(3) Solve the two-point boundary value problem 
using subroutine QUASI. 

(4) If it is desired, determine if the resulting 
solution satisfies the endpoint sufficiency 
condition using subroutine FOCAL. 

(5) Verify the solution determined by QUASI by inte- 
grating the nonlinear equations with the initial 
values found by subroutine QUASI. 

(6) If other solutions are to be found to end- 
points near the solution just determined, use 
the solution just determined as the initial 
solution estimate and return to step 3- If 
not, return to step 1 to solve a new problem. 
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Program SPACE, the first program in Appendix B, follows 
the above outline. 

5-4 Examples of Multiple Stationary Solutions 

In the course of the following discussions, ref- 
erence will often be made to actual orbital transfers 
such as the one shown in Figure 5»4* The interpretation 
of these figures will be aided by the following conven- 
tions : 

(1) The initial orbit is a circular orbit; initial 
motion in this orbit is counter-clockwise. 

(2) Short lines eminating from the transfer trajec- 
tory at selected points as shown in Figure 5*4 
indicate the true direction of the thrust vec- 
tor. The tail of the thrust vector is always 
shown in contact with the transfer trajectory. 

(3) The first thrust vector (proceeding in a counter- 
clockwise direction) on the initial orbit indi- 
cates the point of thrust initiation. The final 
thrust vector indicates the point of thrust 
termination. 

(4) Although the thrust vector is displayed at a 
finite number of points along the transfer tra- 
jectory, the vehicle is thrusting continuously 
from the point of thrust initiation to the point 
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of thrust termination. Likewise the thrust vec- 
tor rotates in a continuous manner. The highly 
nonlinear nature of this rotation causes the 
orientation of the thrust vector to change sud- 
denly at times as shown in Figure 5-^-* 

(5) Motion in the terminal orbit is in a counter- 
clockwise direction, except where noted other- 
wise . 

The computational procedure for implementing the 
endpoint sufficiency test, as described in the previous 
section, can be applied to the multiple stationary solu- 
tions discussed in section 5-2. From Figure 5*3 it can 
be seen that two solutions exist to point B. This point 
can be reached by starting integration with an initial 
control angle of either 4-° or 16 ° . 

Using the modified Newton-Raphson method for 
solving two-point boundary value problems (see Chapter 3)> 
two solutions which terminated at a point near point B 
were found. Both solutions result in a final nondimen- 
sional semimajor axis A of 1.58 and a final eccentricity 
e of .110. The first trajectory begins, with an initial 
angle of approximately l 6 ° ( 15 . 5 °;. This orbital trans- 
fer is pictured in Figure The thrust vector con- 

trol angle increases from 16 ° to over 50 ° and then sud- 
denly decreases rapidly to - 76 °. The control angle then 
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slowly increases to -36°. In this case the state space . 
endpoint is given hy u^, = . 82 6 , g^ = 5 - 92 °, and :R^ = 1 . 52 . 
It Is readily verified that these states correspond to 
A^, = I .58 and e^ = .110. Using the computational proce- 
dure discussed in section 5 - 3 ; the above endpoints were 
found to satisfy the endpoint sufficiency test. 

The second trajectory Initiated transfer with 
x = 4° (4.4°). This orbital transfer is shown in Figure 
5 -5- In this case the thrust-vector control angle 
increases slowly to approximately 80 ° and then suddenly 
swings past l80° to approximately -100°. The control 
angle then increases at a slower rate until x^ = -31° • 

Note that the thrust vector has actually been In a 
"retro-fire" position for a short period of time during 
this maneuver. The trajectory terminates with the fol- 
lowing final status: u = . 773 , g f = 6 . 10 °, and R, = 

1.62. These states , as before, correspond to = I .58 
and e f = .110. Numerical results show that this trajec- 
tory and its resultant endpoints do not satisfy the endpoint 
sufficiency conditions. The corresponding trajectory is 
therefore nonoptimal. In fact, the above nominal end- 
points represent a local maximum value for the func- 


tional J* . 
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This fact was substantiated by computing neigh- 
boring stationary trajectories to fixed final points in 
the vicinity of the nominal final points. These final 
points were chosen beforehand to correspond to A f = 1.5-8 
and e^ = .110. In all cases these neighboring trajec- 
tories had performance indices less than those of the 
nominal trajectory. 

Since T = T from T = 0 until t = the 
max f 

terminal time is proportional to the magnitude of the 

performance index J of equation (5-1.25) • The value of 

is therefore a convenient measure of the performance 

of a given maneuver. For the nonoptimal transfer 

initiated with x = 4°, the value of t„ was found to be 

o f 

5 . 89 . For the optimal transfer initiated with x Q = l6°, 
the value of x- was 4.38. The optimal trajectory there- 
fore yields a performance index which is 25 $ less than 
that for the nonoptimal trajectory. 

Figure 5*6 shows the trajectories in A-e space 
corresponding to the optimal and nonoptimal orbital 
transfers shown in Figures 5-4 and 5*5; respectively. 

A trajectory in A-e space should not be confused with 
the locus of optimal endpoints previously discussed. A 
trajectory in A-e space depicts the values of A and e 
which were attained during the progress of the transfer. 
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ectory to A — 1.58 and e = . 11 
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The loop in the nonoptimal trajectory corre- 
sponds to the sudden rotation of thrust vector to the 
retro-fire position as depicted in Figure 5-5- 

If the thrust magnitude were an allowable con- 
trol, this transfer would obviously be nonoptimal since 
the thrust could then be turned off at point B in Figure 
5-5 and reignited at a later time. Doing so would save 
all the effort needlessly expended traversing the loop. 

It is interesting to note that the argument of periapsis 
for the nonoptimal final orbit (337°) is almost diamet- 
rically opposite to that of optimal final orbit ( 119 °), 
where periapsis is measured in a counter-clockwise 
direction from the point of thrust initiation. 

In terms of the locus of "optimal" endpoints 
shown in Figure 5-2, one endpoint on branch CD has been 
determined to be nonoptimal. The application of the 
endpoint sufficiency test to other endpoints on branch 
CD has shown these endpoints to be nonoptimal, whereas, 
other endpoints on endpoint locus OABC have consistently 
passed the endpoint sufficiency test. The test is inde- 
terminate at point C, since this point represents a 

2 

final circular orbit and sin g^ = 1 - u^R^ = 0 in equa- 
tion (5*3-6). 

Since similar results have been observed for 
loci generated with initial control angles other than 4°, 
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the above conclusions can be generalized to eliminate 
from consideration branches corresponding to branch CD. 
This generalization applies only if the locus of optimal 
endpoints exhibits the same qualitative characteristics 
as those shown in Figure 5*2, With such branches elim- 
inated from consideration, the multiplicity of solutions 
shown in Figure 5-3 disappears. 

A second example of multiple stationary solu- 
tions is illustrated in Figures 5*7 - 5-9- l n this case 
the optimal trajectory shown in Figure 5-7 starts with 
x q = 3-53° and has a thrusting period of 6.10 nondimen- 
sional time units. The nonoptimal trajectory shown in 
Figure 5-8 initiates thrusting with x = - 4.26° and has 
a thrusting period of 7-50 nondimensional time units. 

The latter trajectory was numerically determined to be 
nonoptimal through the use of the endpoint sufficiency 
condition. In this case the optimal trajectory has a 
performance index which is more than 20$ less than that 
for nonoptimal transfer. The qualitative remarks per- 
taining to the first example also apply in this case. 

The main difference between these two examples is a com- 
putational one. Both trajectories in the first example 
were terminated when the first zero of the cutoff func- 
tion was encountered. However, the nonoptimal trajectory 






Fig. 5-9 An Optimal and. a Nonoptimal Trajectory to A = 2.0 and 
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of the second example shown in Figure 5-8 was terminated 
when the second zero of the cutoff function was encoun- 
tered . 

The endpoint sufficiency test is an effective 
computational tool because it allows the rejection of 
certain nonoptimal solutions. When multiple solutions 
are encountered, the true optimal is easily distinguished 
via a direct comparison of the performance indices of 
the two solutions. Using the endpoint sufficiency test, 
a complete class of nonoptimal solutions can be discarded 
immediately upon encounter. Without the aid of the end- 
point sufficiency test an investigator has no indication 
that solutions he is generating are nonoptimal until he 
encounters multiple solutions. Hence, an overall sav- 
ings in computation time is realized by using the end- 
point sufficiency test since the amount of time wasted 
generating nonoptimal solutions is greatly reduced. 

As shown in section 1.3, multiple solutions can 
also arise if path sufficiency conditions are not met. 

As an example. Figures 5 . 10 - 5-12 summarize a case in 
which two orbital transfers exist to A = .835 and 
e = .270, both of which satisfy the endpoint sufficiency 
condition. The optimal transfer initiates thrusting 
with x Q = 47°, but eventually finishes the transfer with 
a period in which the thrust vector is in a retro-firing 
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Fig. 5.12 Two Orbital Trajectories to A = .835 and e = .270 which both 
Satisfy the Endpoint Sufficiency Condition 
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position. The nonoptimal transfer is initiated with a 
retro-firing period followed by a period with a compon- 
ent of the thrust in the forward direction. The anti- 
symmetry of these two transfers is apparent from Figure 
5 - 12 . In this case the optimal trajectory has been dis- 
tinguished from the nonoptimal trajectory by a direct 
comparison of performance indices. The optimal transfer 
is achieved with a thrust duration of 2.71 nondimensional 
time units, while the nonoptimal transfer requires 2.75 
time units. Future investigations can explore the possi- 
bility of eliminating multiple solutions of this type in 
a direct manner by implementing a path sufficiency test. 

5«5 Optimal Orbital Transfer Results 

This section presents further results concerning 

the minimum fuel orbital transfer problem described in 

section 5*1* The investigation was again restricted to 

problems in which the thrust magnitude was constant. In 

addition only a single set of space vehicle parameters 

were investigated. A nondimensional thrust F = .06666667 

+50 

and an exhaust velocity u g = 10 ^ were used. The thrust 
was chosen to coincide with that used by Mason (1967)* 
Rather than choose a single arbitrary value for the 
exhaust velocity, it was set equal to the limiting value 
of infinity. Using the modified Newton- Raphs on method 
to solve the two-point boundary value problem (see Chapter 
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3) and the endpoint sufficiency test to aid in the deter- 
mination of the true optimality of solutions, optimal 
trajectories were computed to orbits having a wide range 
of final semi-major axis and eccentricities. 

Figure 5-13 shows the loci of optimal endpoints 
determined for various initial thrust-vector control 
angles x . There are no intersections among these loci. 
As mentioned earlier, various regions of this endpoint 
space are attained by terminating the thrust at various 
zeros of the cutoff function. Figure 5»l4 depicts 
regions in which the first, second, third, and fourth 
zero have been used. Crosshatched areas indicate regions 
of uncertainty. Final orbits represented by points in 
the region labeled "l" are attained by terminating thrust 
when the first zero of the cutoff function is encoun- 
tered. Final orbits represented by points in regions 
labeled ,f 2" are attained by terminating thrusting after 
the second zero of the cutoff function has been encoun- 
tered. Similar statements can be made concerning the 
regions labeled "3" and "4. 1 ' 

Some of these boundaries seem to have little 
physical significance. For example, a comparison between 
trajectories in the regions labeled A and those in the 
region labeled B shows there to be no sudden discontin- 
uity of any kind across the boundary. The same is true 
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Fig. 5.14 Optimal Number of Cutoffs to Attain Various Orbit 
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for the boundary between regions C and D. This is not 
true, however, for the boundary between regions B and C. 
This fact is apparent from Figure 5-13- In this figure 
the boundary between regions B and C lies just above the 
line with positive slope labeled -26°. From Figure 5-13 
it is seen that there is a jump in initial angle x Q across 
this boundary, viz., -32° to -2 6°. 

An explanation for this jump is shown in Figure 
5.15* Curve A shows the cutoff function for a trajec- 
tory initiated with x Q = -32° and x Q = .46. Thrusting 
is terminated at the second zero of this function. This 
corresponds to an initial orbital transfer to a point at 
the end of the curve labeled -32 in Figure 5-13- Curve 
B in Figure 5-15 shows the cutoff function for a trajec- 
tory initiated with x Q = -32° and x Q = .44. Although 
only a slight change has been made in the cutoff func- 
tion, the second cutoff does not occur, and the corre- 
sponding orbital transfer is to a hyperbolic orbit. 

This accounts for the physical discontinuity of orbital 
transfers to points just on either side of the boundary. 

Figure 5-16 depicts the transfer time needed to 
attain an entire range of various final orbits . The 
transfer time is proportional to the performance index 
of equation ( 5 - 1 . 25 ) • To attain a given final serai- 
major axis, this figure shows it is easier to reach an 
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orbit with an intermediate eccentricity than it is to 
reach a circular orbit or an orbit with high eccen- 
tricity. The dashed line in Figure 3.16 represents the 
locus of final eccentricities associated with minimum 
effort transfers to any given semi-major axis provided 
no constraint is placed on eccentricity . 

Two fundamental principles of astrodynamics can 
be verified by observing Figure 5-17- This figure pic- 
tures an orbital transfer to a final orbit with a semi- 
major axis of 1.5 and an eccentricity of . 60. This is a 
typical orbital transfer requiring a substantial change 
in both energy and angular momentum. From the figure 
it can be seen that the thrust vector is essentially 
aligned with the velocity vector during the portion of 
transfer nearest the center of attraction. It is well 
known that a small change in the vehicle's kinetic energy 
is given by 

AE = mv-Av^ (5.5-1) 

Therefore , the energy of the vehicle is increased in the 
most economic fashion by thrusting in the direction of 
the vehicle velocity when the vehicle velocity is a maxi- 
mum. For low thrust orbits this generally occurs at 
points on the transfer trajectory nearest to the center 
attraction. In Figure 5-17 the space vehicle is using 
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its thrusting capability to best advantage by increasing 
its energy during the initial portion of the orbit . 

On the other hand a change in angular momentum 
is most economically accomplished by thrust perpendicu- 
lar to the radius vector at the local apogee of the 
transfer. This is easily seen from the following equa- 
tion 

Ah = mr X Av 

The product r X Av is maximized by thrusting in a direc- 
tion perpendicular to the radius vector at points along 
the transfer trajectory at which r is a maximum. 

It will be seen later that an optimal maneuver 
for making large changes in eccentricity , but not in 
energy, requires that the vehicle first make a large 
gain in energy. Then at large radii the eccentricity is 
economically changed. Figure 5.18 shows an orbital 
transfer yielding an extreme change in angular momentum. 
In fact, the sense of the angular momentum vector 
reverses during the transfer and orbital motion changes 
from counter-clockwise to clockwise. 

Figures 5-19 ~ 5*21 illustrate several optimal 
orbital transfers to final circular orbits. The trans- 
fer shown in Figure 5-19 is typical of optimal transfers 
to final circular 'orbits with 1.0 > R > I. 98 . These 
transfers are characteriz.ed by a sudden rotation of the 
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thrust vector through a retro-fire position. This often 


generates loops in the A-e space trajectories (see 
Figure 5-21) which have previously been associated with 
nonoptimal transfers. As previously mentioned,, the end- 
point sufficiency condition is indeterminate for trans- 
fers to circular orbits. Loops in A-e trajectories may 
be a result of requiring the thrust magnitude to be con- 
stant. For near orbit transfers a Homan impulsive trans- 
fer (Lawden, 1963, pp. 106-110) proves to be more optimal 
and one may conclude that thrusting at intermediate radii 
is to be avoided under certain circumstances. In a sense 
the loop represents a period in which the thrust should 
be off. 

The transfer shown in Figure 5*20 is typical of 
optimal transfers to circular orbits with R > 2.00. In 
these cases the thrust vector exhibits a highly nonlinear 
oscillation about x = 0„ but never passes through the 
retro-fire position. In Figure 5-21 notice that trajec- 
tories to circular orbits first pick up a substantial 
portion of the final energy required by proceeding along 
the path of minimum effort. The eccentricity is then 
adjusted to zero by proceeding almost perpendicular to 
lines of constant effort. 

Figures 5-20, 5-22 and 5-23 summarize transfers 
to final orbits with a semi-ma.jor axis of 3-0. The 
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fundamental principles of astrodynamics are especially 
evident in Figure 5-22. Energy is first increased near 
the center of attraction by aligning the thrust vector 
near the velocity vector. Eccentricity is then increased 
by thrusting nearly perpendicular to the radius vector 
at large values of the radius . 

Figure 5-24 shows several optimal orbital trajec- 
tories to a final eccentricity of . 90. All of these 
transfers are initiated by a large gain in energy attained 
by proceeding along the minimum effort path. The space 
vehicle then increases its eccentricity at points far 
from the center of attraction. These trajectories again 
represent excellent examples of the fundamental prin- 
ciples of astrodynamics. 

Figures 5'- 25 - 5-28 illustrate optimal orbital 
transfers to final orbits with a semi-major axis of one. 

In all cases only a change in final angular momentum is 
required. However, in accordance with the second funda- 
mental principle, such a change is most economically 
performed at large radii. Consequently each of these 
transfers is initiated with an effort to increase the 
energy of the orbit. From Figure 5-28 it can be seen 
that this is done by proceeding along the path of minimum 
effort. When large radii are achieved, the thrust vector 
is pointed nearly perpendicular to the radius vector in 





171 



Fig. 5-2 6 


An Optimal Orbital Transfer 
to A = 1.0 and e = . 50 
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an effort to achieve the final angular momentum desired. 
At the same time a component of the thrust vector exists 
in the retro-fire position and the semi-major axis is 
decreased hack to one. 

Figures 5-29 - 5-32 illustrate orbital transfers 
to final orbits with semi-major axes of 0.75* A com- 
pletely different region of transfers is encountered. 

To attain final orbits below the minimum effort curve 
(see Figure 5-32), the thrust vector is initially posi- 
tioned in a retro-fire position and oscillates about 
x = 180°. Two representative optimal transfers are shown 
in Figures 5-29 and 5-30. 

In transferring to inner orbits with high eccen- 
tricity , it is difficult to lose the angular momentum 
required to achieve highly eccentric orbits , since the 
final semi-major axis must be small. In other words , 
the fundamental principle for losing energy is at odds 
with the principle for decreasing angular momentum. To 
achieve orbits with high eccentricity , the change of 
angular momentum required is the dominant factor. Con- 
sequently, to attain final orbits with eccentricities 
above the minimum effort path, the thrust vector is 
initially positioned with a component of thrust in the 
forward direction. The vehicle energy is first increased 
and as large radii are encountered, the energy and 
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angular momentum are simultaneously decreased. Such a 
transfer is shown in Figure 5-31- In this case, thrust 
expended near the perigee of transfer is largely lost . 
If thrust magnitude were a control, the thrust would 
certainly have been turned off during this portion of 
transfer . 



CHAPTER 6 


CONCLUSIONS 

Multiple stationary solutions are frequently 
obtained with optimization problems using the calculus 
of variations. The existence of multiple stationary solu- 
tions can be attributed to one of several distinct rea- 
sons. A reason that multiple solutions are frequently 
obtained is that the complete set of necessary and suffi- 
cient conditions for a solution have not been imposed. 

It has been shown here that a sufficiency test 
for the Problem of Bolza can be broken down into two inde- 
pendent tests: 

(1) a path sufficiency test, and 

(2) endpoint sufficiency test. 

An endpoint sufficiency test is developed here for prob- 
lems with variable endpoints. If a solution candidate 
satisfies both test as well as the first three necessary 
conditions, it simply satisfies a sufficiency condition. 
That this sufficiency condition is also necessary has not 
been proven. 

Multiple solutions can also occur in problems with 
singular control. To date the complete set of necessary 

1 8o 
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and sufficient conditions needed to distinguish the true 
optimum have not been formulated. 

Problems also exist for which the calculus of 
variations provides no criteria for selecting the global 
minimum. The theory of the calculus of variations is 
effective only in determining trajectories which are local 
optimums. Only small variation in the state space trajec- 
tory are considered. Problems with bounded control and 
periodic solution often have several local optima. 

It has been shown that once the first necessary 
path conditions have been applied , a calculus of varia- 
tions problem with variable endpoints is reduced to a 
problem of the minimization of a function of several 
variables . Analytical application of the endpoint suf- 
ficiency condition requires the analytical integration 
of the set of state variable and Euler-Lagrange differen- 
tial equations. Since in most cases this is difficult 
or impossible, an algorithm has been developed for the 
numerical Implementation of the endpoint sufficiency 
test . 

The endpoint sufficiency condition has been 
shown to be an effective computational tool In complex 
applications. For example, through the use of the end- 
point sufficiency test, a complete class of nonoptimal 
solutions can be discarded immediately upon encounter. 
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Without the aid of the endpoint sufficiency test an 
investigator would have no indication that solutions he 
is generating are nonoptimal until he encounters multiple 
solutions. This results in a. great savings in computa- 
tion time. 

The alternative of determining the set of all 
multiple stationary solutions to a problem, so that the 
true optimal may be distinguished via a direct compari- 
son of the performance indices, is often a formidable 
task. No criteria exist for determining the number of 
multiple solutions that exist and an investigator could 
never be sure he had found all of the multiple solu- 
tions . 

The difficulty of solving a two-point boundary 
value problem needed to Implement the endpoint sufficiency 
condition was eliminated by using the generalized Newton- 
Raphson algorithm. It has proved to be a powerful com- 
putational tool. 

The comprehensive optimal orbital transfer 
example gave some interesting insights Into the general 
problem of low thrust minimum fuel transfers . Qualita- 
tive aspects of all orbital transfers were found to be 
consistent with fundamental principles postulated for 
changing energy and angular momentum. 



APPENDIX A 

MINIMIZATION OF A FUNCTION OF SEVERAL VARIABLES 

This appendix presents a rigorous proof of the 
necessary and sufficient conditions for a function of 
several variables to be a minimum when subject to alge- 
braic equations of constraint. These proofs could be 
established easily, but with less rigor , through the use 
of Lagrange multipliers as was done in section 1.2. The 
validity of using such multipliers in establishing neces- 
sary conditions has been rigorously verified. However , 
the validity of using Lagrange multipliers to establish 
sufficiency conditions in control notation has not been 
established with rigor until recently (Vincent, 1969) • 

In the proof that follows, free use will be made 
of the notation and conventions established at the begin- 
ning of section 2.3. Following the methods of Vincent 
(1969), consider the problem of minimizing a function of 
several variables 


J = J(w,v) (A. 1) 

subject to the constraints 

A (w,v) ~ A (A* 2 ) 
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where both J and. jr_ are functions of class C and the con- 
straints are such that the determinent of the Jacobian 


8v 


(A.3) 


is nonsingular. The dimension of Jj_ and w is assumed to 
be p and the dimension of v is assumed to be q. 


A.l Method of Implicit Functions 

Since is of C 2 and condition (A.3) has been 
postulated, the implicit function theorem (Buck, 1965, 
pp. 283-286) states that equation (A. 2) implicitly 
assures the existence of the vector function ¥ explic- 
itly relating the dependent variables w to the independent 
variables v 


^(v) 

Wg(x) 

Wp(x) 


(A.4) 


By substituting (A. 4) into (A.l), J becomes a function 
of v only 

J = J(¥(v),v) (A. 5) 

Define the general value of independent variables 
v in a small neighborhood of an optional point v° 



V 


o , 

= v + e c 
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(A. 6) 


where c is a vector of' arbitrarily chosen small, but non- 
zero^ constants and e is a scaler multiplier. Then, from 
(A. 1) and (A. 2) 


J = J [¥(v° + e _c), v° + e cj (A. 7) 

i [¥(v° + e _c ) , v° + e cj = 0 (A. 8) 


Now J is a function of e only, and the necessary 
condition for an ordinary local extremum is 


dJ 

de 


Sj T h 


rp 

+ £ 
chF 


= 0 


(A. 9) 


d¥. 

where h. is the vector with elements h. = — . The 

vector h represents changes in the dependent variables 
w corresponding to the changes c_ in the independent vari- 
ables . Differentiating equation (A. 8) with respect to e 
yields 


' chjj_‘ 
^w 


h -i 



0 


Solving for h yields 



-1 

- djr 

ciW 


3V 


(A. 10) 


(A. 11) 


Substituting (A. 11) into (A. 9) and rearranging gives 
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dJ _ 

8 J T 8j T 


-1 

‘ 8^- 


de 

37 ~ 

Sw 


_37 



(A. 12) 


Since c_ is an arbitrary nonzero vector, each of the ele- 
ments of the vector in parenthesis must be equal to _0 at 
the optimal point . 


8J 

3v 


T 


8j 

ciw 


T 

' a±“ 

-i 

~ 8 ^ ' 


37 


3v 



(A. 13) 


If equation (A. 13) is satisfied, then a sufficient 
condition for a 'local minimum is that 


dfj 

de 2 v 


(a.i4) 


o 


must be positive definite for arbitrary values of h and 
c_ satisfying equation (A. 11). 

Before evaluating this expression, an identity 
for taking the partial derivative of the inverse of a 
matrix must be developed. Let A . . represent a general 


element of the 


'5w 


-1 


matrix : 


A ij 


^1/ 

"5w 


-1 

ij 


(A. 15) 


Then in indicial notation the definition of 
inverse may be expressed as 

6 . = A . 

qj % mj 


(A. 16) 
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where 8 . is the Kroniker delta. Premultiplying by A. 

0.0 10 


gives 


A ij A iq ^0 


(A. 17) 


Taking the partial derivative of both sides yields 


8 


d 


d~r„ ^ A i j ^ dr^ ( A iq^ + A iq ’Sr’ 


' n 


' n 


n 


q 


A . 


i^T/ mj 


+ ^ im dr„ ^ A m j ^ 


(A. 18) 


n 


Since 8. . represents a constant, equation (A.l8) reduces 
10 


to 


d 


dr n ^ A i j ^ ‘ "5r n ^ A ij ^ + A iq ’dr' n 


A . - 3 — — 

ov mj dr^ v ”i0 


(A,,) 


i m 


(A. 19) 


which gives the desired identity: 


or 

n 


(A if 




A . 

mj 


(A. 20) 


Prom this point on results must be expressed in 

^ I dijf \ 

indicial notation since - 5 — -y-2. is a tensor. In 

dr dw 
n \ ml 

indicial notation equation (A. 12) becomes 


dJ 

de 


5v^ ~ 3w7 A ij dv k j 


dJ dJ 


(A. 21) 



Using equations (A. 20) and (A. 21) the sufficiency con- 
dition (A. 14) yields 
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O . p p 

d J d^J h d^J 

de 2 9 V T k k 1 d V v k k n 


0 2 J -^ili •; -i 

- q y a iL ph - A - c c 

3w_Sw. « kn dv n dw i k n 


n i 


~n t d ■ diji • 

JL a ... 9 , A . J. q 

8w i iq dw n dw~ i rnj dv^ k n 


(A. 22) 


dj 


dw. A iq dv dw ^raj dv, u k u n 


d-llf 


a 


dilf . 

A . -s-~ c, c 


n m 


k 


dj 


dw7 A ij dw p 5v' k c k n p 


d^ , 


c, h 


dJ n 

dw7 A ij dv n dv k C k C n 


must be positive definite for arbitrary values of 
£ and h satisfying equation (A. 11). Recognizing the 
indicial notation representation of equation (A. 11) 

d ijf 1 

h m = d V “ C k (A. 23) 

in four terms of equation (A. 22) and regrouping terms , 
condition (A. 23) can be written 
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' w v i 

n k 


A 

chfT A ij 


11 
r 6v, 
n k 


c, c 
k n 


8 2 J_ 8j 
"5w^5v k “ pj 


y 

r . 5vT 
1 k 


r dw.' dw pq dv dw. 
n i p ^ n x 


U U U U « 

8m .dw. dw pq 
Ji P 


3w . dw . 
i 0 


c, h. 
k i 


h. c 
l n 


h. h . 

i 0 


(A. 24) 


must be positive definite for arbitrary values of c 
and h satisfying equation (A. 11). 

Equation (A. 24) and equation (A. 13) provide a 
set of necessary and sufficient conditions for J (w, v) to 
be minimum when the equations of constraint jr(w,vj = 0 


must be satisfied, 


A. 2 Method of Lagrange Multipliers 

The necessary and sufficient conditions can be 
put in a form which is more convenient to use by defining 
the augmented function 

J*(w,v,^j - J(w,v) + ^ f(w,v) (A. 25) 

where ja is a vector of constant multipliers called 
Lagrange multipliers . 

If the Lagrange multipliers are given the identity 

dJ A 

^ j <5wT A ij 


(A. 26) 
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or in matrix notation., 


J± 


dj T ' 


-i-i 




(A. 27) 


several observations can be made. Thus necessary condi 
tions for J to be a minimum (A. 13) become 

dj* 


dv 


0 


(A. 28) 


The sufficiency conditions (A. 24) may be written as 


/ 




d 2 J 




j - T ^ 

dv n dv k. + ^ j <3v n 8v~ ) c k c n + I ^w73v^ + ^ j '5w75v 


i k 


+ 


1 b 2 j a *4 

d V w k ^ 6v n dw i 


d 2 J a2 f q ' 

^i c n 4 1 cihldhT + ^q dw.. dw 

J J- 


c k h i 

(A. 29) 

\ 


h. h . 
i J 


i J j 

must be positive definite 

for arbitrary c and h satisfying equation (A. 11). Define 
the vectors r and d as 


W 1 


h i 

w 2 


h 2 

w p 

d - 

*P 

V 1 


C 1 

v 2 


C 2 

V 

q 


c 

q 


(A. 30) 
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The sufficiency condition (A. 29) may now be put 
in compact matrix notation by using equation (A. 25) and 
definitions (A. 30): 


d 


T 


d 2 J* 

chFSr 


(A. 31) 


must be positive definite 

for arbitrary _c and h satisfying equation (A. 11). To 
guarantee that c. and h in vector _d satisfy equation 
(A. 11), elements of d. may be expressed as functions of 
the independent constants c only. In matrix notation 
this may be expressed as 


where 



‘ djT 

-1 

" cH 

_3w 


3V 


(A. 32) 


(A. 33) 


and I is a q by q identity matrix. The sufficiency con- 
dition (A. 32) may then be written as 


T T 
c Q 


d 2 J*' 

drdr 


Q c 


(A. 34) 


must be positive definite 

for arbitrary values of the elements of c . 

The advantage of the Lagrange multiplier tech- 
nique is that the necessary and sufficient comditions can 
be expressed in compact and easy-to-remember matrix nota- 


tion. 



APPENDIX B 
PROGRAMS 
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APPENDIX B -A 1 


A 2 
A 3 

4 4 

* ORBITAL TRAJECTORIES * A 5 

❖ A 6 

A 7 

MAINLINE PROGRAM A 8 

PROGRAM SPACE CONTROLS THE GENERATION OF SOLUTIONS TO AN OPTIMAL A 9 

PROBLEM IN CELESTIAL NAVIGATION. THE PROBLEM IS TO TRANSFER A LOW A 10 

THRUST SPACE VEHICLE FROM AN INITIAL CIRCULAR ORBIT TO ANY COPLANAR A 11 

ELLIPTIC ORBIT OF GIVEN ENERGY AND ANGULAR MOMENTUM. THE THRUST A 12 

PROGRAM IS TO BE DETERMINED SO AS TO MINIMIZE THE INTEGRAL OF THE A 13 

PRODUCT OF THRUST MAGNITUDE AND TIME. THE FINAL ARG. OF PERIAPSIS A 14 

THE FINAL TRUE ANCMOLY, THE FINAL RANGE ANGLE, AND THE FINAL TIME A 15 

ALL ARE UNSPECIFIED AND CONSIDERED FREE. A 16 

A 17 

THIS PROGRAM AND ASSOCIATED SUBROUTINES ARE EXPLAINED IN DETAIL IN A 18 

CHAPTERS 5 AND 3 AND IN SECTION 2.7 OF THIS DISSERTATION. A 19 

INPUT A 20 

A 21 

EACH SET OF INPUT DATA- REQUIRES TWO CARDS, A 22 

CARD 1 A 23 

XI = ESTIMATE OF INITIAL THRUST CONTROL ANGLE. A 24 

XDOTI = ESTIMATE OF INITIAL NOND I MEN S I ON AL THRUST CONTROL A 25 

ANGLE RATE. A 26 

LRI = ESTIMATE OF INITIAL RADIUS LAGRANGE MULTIPLIER. A 27 

NEEDED ONLY FOR PROBLEMS WITH MULTIPLE THRUSTING ARCS. A 28 

AF = FINAL NOND I MENS I ON AL SEMIMAJOR AXIS DESIRED. A 29 

NEEDED FOR FOR BOUNDARY VALUE PROBLEMS ONLY A 30 

= 0 IF THE FLOODING TECHNIQUE IS BEING USED AND NO A 31 

BOUNDARY VALUE PROBLEM is TO BE SOLVED. (SECTION 3.1.1) A 32 

ECCF = FINAL ECCENTRICITY DESIRED (FOR BOUNDARY VALUE PROBLEM) A 33 

TOL = MAXIMUM INTEGRATION ERROR/ UNIT STEP. (APPROXIMATE ONLY) A 34 

CTO = A CONVERGENCE TOLERANCE FOR TERMINATING ITERATION IN A 35 

SUBROUTINE QUASI CURING SOLUTION OF A BOUNDARY VALUE PR. A 36 

PARDEL = PERCENTAGE CHANGE IN FINAL ENDPOINTS WHEN NUMERICALLY A 37 

COMPUTING PARTIAL DERIVATIVES FOR ENDPOINT SUFFICIENCY A 38 

TEST. (SEE SECTIONS 2.7 AND 3.3) A 39 

A 40 

CARD2 A 41 

NS = RADIAL VELOCITY DIRECTION FOR ELLIPTICAL INITIAL ORBITS. A 42 

= -1 IF APROACHING PERIAPSIS A 43 

KFREQ = MAXIMUM NUMBER OF INTEGRATION STEPS BETWEEN PRINTED A 44 

TRAJECTORY OUTPUTS. A 45 

NV = NUMBER OF THE VARIABLE WHICH IS TO BE CHANGED BEFORE A 46 

AUTOMATICALLY REPEATING A SOLUTION. A 47 

= 0 IF A NEW SET OF DATA CARDS IS TO BE READ IMMEDIATELY A 48 

FOLLOWING THE SOLUTION GENERATED BY THIS CARD. A 49 

= 1 AND AF= 0, CHANGE XI BY DV. A 50 

= 2 AND AF= 0, CHANGE XDOTI BY DV. A 51 

= 3 AND AF= 0, CHANGE LRI BY OV. A 52 

= 4 AND AF= 0, CHANGE FI BY DV. A 53 

= 1 AND AF NOT EQUAL TO 0, CHANGE AF BY DV. A 54 

= 2 AND AF NOT EQUAL TO 0, CHANGE ECCF 8Y DV . A 55 

= 3 AND AF NOT EQUAL TO 0, CHANGE FI BY DV. A 56 

= 4 AND AF NOT EQUAL TO 0, CHANGE UE 8Y DV. A 57 

NCH = NUMBER OF CHANGES, DV, TO BE MADE IN VARIABLE NV ABOVE. A 58 

NT YP =1 IF AN INITIAL CIRCULAR ORBIT AND A SINGLE THRUSTING A 59 

PERIOD ARE DESIRED. A 60 

= 2 IF AN INITIAL CIRCULAR ORBIT AND MULTIPLE THRUSTING A 61 

PERIODS ARE DESIRED. NOT fully OESU&oeD- A 62 

= 3 IF AN INITIAL ELLIPTIC ORBIT WITH MULTIPLE THRUSTING A 63 
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PERIODS IS DESIRED. THIS OPTION IS MOT FULLY DEBUGGED. A 64 

IF = I IF AN ENDPOINT SUFFICIENCY CONDITION TEST IS TO BE A 65 

PERFORMED. A 66 

NB = NUMBER CF CUTOFFS. THE SOLUTION IS TERMINATED WHEN THE A 6,7 

NB-TH ZERO OF THE CUTOFF FUNCTION IS ENCOUNTERED. A 68 

ID = DUMMY FOR FUTURE EXPANSION. A 69 

NIT = MAXIMUM NUMBER OF GENERALIZED NEWTON-R APHSON ITERATIONS A 70 

ALLOWED FOR CONVERGENCE OF A SOLUTION TO THE BOUNDRY A 71 

VALUE PROBLEM. A 72 

DV = VALUE OF THE CHANGE TO BE MADE IN VARIABLE -NV- ABOVE. A 73 

PCT = A MAXIMUM PERCENTAGE CHANGE ALLOWED IN INITIAL VALUES A 74 

DURING THE SOLUTION OF A BOUNDARY VALUE PROBLEM. A 75 

THIS PARAMETER IS USED TO CONTROL THE STABILITY OF THE A 76 

CONVERGENCE OF THE SOLUTION. (SEE SUBROUTINE QUASI! A 77 

FI = MAXIMUM NONDIMENSICNAL THRUST MAGNITUDE. A 78 

UE = NONDIME NSI ON AL EXHAUST VELOCITY. A 79 

A 80 
A 81 

PROGRAM SPACE! INPUT, OUTPUT, PUNCH, TAPE 2=1 N PUT, TAPE3=0U TP UT,TAPE4=PU A 82 

INCH) A 83 

REAL Y! 100) , P( 20 ) ,CONST( 10) , X0LDI6, 1000) , YINITt 10! A 84 

REAL M,LU,LG,LR,LM, LRI A 85 

COMMON /MAIN/ E , H, S MA , ECC , CUT , S , APO , PER , X , W I ER , DPR , HA M A 86 

COMMON /INTV AR/ T I ME , Y , P , NE , TOL , DX , MODE A 87 

COMMON / 1 OLD/ XOLD A 88 

EQUIVALENCE (Y(1),U), (Y(2),G), (Y(3),R), (Y(4),M), (Y(5),LU), (Y{ A 89 

16), LG), ( Y ( 7 ) , LR ) , (Y18),LM), (Y(9),AN), lY(iO).TE) A 90 

EXTERNAL DIFFI, BOUND, NOMLIN A 91 

EXTERNAL F I X E D , CONST RT A 92 

EXTERNAL S W I TCH , CUTOFF A 93 

DATA DPR/57.2957795131/ A 94 

DATA PI/3.141592653589793/ A 95 

A 96 

READ FIRST INPUT CARD. SEE COMMENTS ABOVE. A 97 

READ (2,2! X I , XDOT I , LR I , AF , ECCF , TOL , CTO , P ARDEL A 98 

FORMAT (8E10.4) A 99 

A 100 

READ SECOND INPUT CARD. SEE COMMENTS ABOVE. A 101 

READ (2,3) NS ,NC , KFREQ , NV.NCH, NTYP, i F ,N8 , I D» TO* NI T , ID ,DV ,PC T , F I , UE A 102 

FORMAT (515,511, 215, 2F5.0,2F10.0) A 103 

IF (FI.LE.O.O) STCP 1111 A 104 

NCS=NC A 105 

N8S=NB A 106 

I D=2 A 107 

A 108 

SET INITIAL DATA ASSUMING A CIRCULAR ORBIT WITH A SINGLE A 109 

THRUSTING PERIOD SPECIFIED. A 110 

U=1.0 A 111 

G=0.0 A 112 

R= 1 . 0 A 113 

M= 1 . 0 A 114 

LM=0.0 A 115 

E=-.5 A 116 

AN=0. 0 A 117 

TE=0. 0 A 118 

SMA=1 . 0 A 119 

X=XI A 120 

XP=X A 121 

XDOT=XDQT I A 122 

LR=LRI A 123 

F=F I A 124 

X=X/DPR A 125 

Y { 1 1 ) =0 . 0 A 126 

Y ( 1 2 ) =0. 0 A 127 

KSW=0 A 128 

UEST=1 . OE +50 A 129 
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FACTOR=~F*SlN(X)/M-U*XDOT+U*U*SiN(X)*SIN(G+X) /R+CO S ( G ) / ( R*R ! A 130 

IF ( NT YP. GT . 3 ) STOP A 131 

GO TO (5,6,7!, NTYP A 132 

A 133 

CIRCULAR INITIAL ORBIT WITH NO MASS RATE AND A SINGLE THRUST. A 134 

X AND XDOT DETERMINE LU, LG, AMD LR. LM IS NOT NEEDED. A 135 

LU=M*C0S(X) A 136 

LG=U*LU*SINIX!/C0S1X) A 137 

LR=LG*FACTOR/IU*U*SIN!X)*COS(G+X)) A 138 

GO TO 8 A 139 

A 1 AO 

CIRCULAR INITIAL ORBIT WITH FINITE MASS RATE. A 141 

X, XDOT AND LR DETERMINE LU, LG, AND LM. A 142 

LG=LR*U*U*SIN!X) *COS(G+X ) /FACTOR A 143 

LU=LG*COS(X)/(U*SIN(X! ) A 144 

LM=(LU*( F*CCS( X)/M)+LG*(F*SIN!X1/(M*U)- 1.0/1 R*R*U ) +U/ R ) - F ) *UE /F A 145 

GO TO 8 A 146 

A 147 

ELLIPTIC INITIAL ORBIT - START INTEGRATION AT PER1APSIS. A 148 

H DETERMINES ORBIT SHAPE SINCE SMA IS NORMALIZED TO 1. A 149 

ECC=SQRT( 1.0-H*H/SMA) A 150 

R=SMA*( 1.0-ECC) A 151 

U=H/R A 152 

GO TO 6 A 153 

P ( 1 ) =F A 154 

P (2 ) =UE A 155 

S=SWITCH(TINE,Y,P) A 156 

IF (StO. 00001.LT. 0.0! PI 1 1 = 0-0 A 157 

CUT=CUTOFF(TIME,Y,P) A 158 

H=R*U*C0S(G! A 159 

A 160 

SAVE INITIAL VALUES FOR USE IN B.V. PROBLEM. A 161 

DO 9 1=1,10 A 162 

YINITf I ! = Y ( I ) A 163 

A 164 

PRINT OUT PROBLEM PARAMETERS. A 165 

XD=X*DPR A 166 

WRITE (3,10) XD, XDOT, LR , F , H , UE , TOL , NS, NC, NTYP, NB A 167 

0 FORMAT (1HI,40X,21H***NEW TRAJECTORY*** //15H X, XDOT, LR = F13.8, A 168 

12F 13.8, 30H THRUST, MOM, EXH VEL, TOL =4F9 . 4, 2X , 4 1 2 1 A 169 

WRITE (3,11) A 170 

1 FORMAT ( Z/52H X/ANGLE ENERGY/A MOM/ECC VEL/LU G/LG 5 A 171 

1 OHRAD/LRMASS/LMSWCH/ CUTH/W I ERA PQ/ PER// ) A 172 

CALL OUTPUT A 173 

IF (AF.NE.O) GO TO 12 A 174 

PUNCH 36, F,UE,CNO,X,X DOT ,U,G»R,LU,LG,LR,LM A 175 

CNO=NC A 176 

A 177 

INITIALIZE 800KEEP I NG PARAMETERS FOR INTEGRATION. A 178 

12 KE RR=0 A 179 

OS=0.0 A 180 

0CUT=0. 0 A 181 

KF=KFREQ A 182 

NB=N8S A 183 

NC=NCS A 184 

T 1 ME=0. 0 A 185 

NE= 12 A 186 

M0DE=0 A 187 

XNB=NB A 188 

NLOOP=10.*XN6/TOL A 189 
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A 190 

INTEGRATION LOOP WITH TESTS FOR END OF TRAJECTORY, THRUST A 191 

CONTROL, PHYSICALLY IMPOSSIBLE SOLUTIONS, AND OUTPUT CONTROL. A L92 
START INTEGRATION LOOP £#£###**❖**#*##<=## A 193 
DO 23 MQ= 1 , NLCOP A 194 

CALL RUNKUT ( T I ME , Y , P, NONL I N , NE , TOL , DX , MODE ) A 195 

A 196 

CHECK FOR SIGN CHANGE IN CUTOFF. IF A CHANGE HAS OCCURRED, A 197 

CONVERGE STATE VARIABLES TO EXACT CUT OFF VALUE. A 198 

CUT=CUTOFF!TIME,Y,P) a 199 

IF (0CUT.EQ.0.0.0R.CUTS0CUT.GT.0.0.0R.MQ.LE.3) GO TO 15 A 200 

CALL CONVERG ( CUTOFF, OCUT, CUT , DX, I .OE-08, NONL I N > A 201 

M0DE=0 A 202 

IF (CUT. EQ. 0.0) 13,16 A 203 

A 2 04 

DETERMINE IF TRAJECTORY IS TO BE ENDED AT THIS ENDPOINT. A 205 

3 NB=NB~1 A 206 

WRITE (3,11) A 207 

'CALL OUTPUT A 208 

WRITE (3,24) A 209 

IF (KS W. EQ. 1. AND. ABS(U -UESTl.LT. 0.02) GO TO 14 A 210 

IF (KSW.EQ.l) NB=NB+1 A 211 

IF (KSW.EQ.O.AND.NB.LE.O) 14,15 A 212 

4 KE RR= l A 213 

GO TO 25 A 214 

A 215 

CHECK FCR SIGN CHANGE IN SWITCH. IF A CHANGE HAS OCCURRED, A 216 

CONVERGE STATE VARIABLES TO EXACT SWITCHING VALUE. A 217 

5 IF (NC.LT.O) GO TC 20 A 218 

S=SWITCH(TIME,Y,P> A 219 

IF (OS.EQ.O.O.OR.S*OS.GT.O.O.OR.MQ.LE.3) GO TO 20 A 220 

SNEW=S A 221 

CALL CCNVERG ( SW ITCH , OS, S , DX , 1 .0E-08 , NONL I N ) A 222 

MODE=0 A 223 

IF (S.EQ.O.O) 17,16 A 224 

6 KEUR=2 A 225 

GO TO 25 A 226 

A 227 

TURN THRUST ON OR GFF DEPENDING ON SWITCH GOING + OR -. A 228 

7 NC=NC-1 A 229 

IF (NC.LT.O) GO TC 14 A 230 

IF ( SNEW.GT . 0.0 ) 19,18 A 231 

8 . P ( 1 1=0.0 A 232 

WRITE (3,11) A 233 

CALL OUTPUT A 234 

GO TO 20 A 235 

9 P ( 1 ) =F A 236 

WRITE (3,11) A 237 

CALL OUTPUT A 238 

A 239 

OTHER HALTS BECAUSE OF PHYSICALLY IMPOSSIBLE SOLUTION. A 240 

0 IF (M*E.LT.O.O) GO TO 21 A 241 

KERR=4 A 242 

GO TO 25 A 243 

A 244 

PRINT OUT AFTER EVERY *KF* INCREMENTS OR AFTER EVERY TEN A 245 

DEGREES CHANGE IN X, WHICHEVER COMES FIRST. A 246 

1 X=ATAN(LG/(U*LU) ) A 247 

KF=KF-1 A 248 

IF (ABS(X-XP).LT. .174533. AND. KF.GT.O) GO TO 22 A 249 

KF = KFRE Q A 250 

XP=X A 251 

CALL OUTPUT A 252 

C A 253 

22 OCUT=CUT A 254 

OS=S A 255 

23 CONTINUE A 256 

KERR=1 0 A 257 

C A END INTEGRATION LOOP A 258 
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C A 259 

C FINAL SUMMARY OUTPUT AFTER EACH SOLUTION. A 260 

24 FORMAT ( 1H ) A 261 

25 FTA=ACOSt (SMA#t 1 . 0-ECC**2 > /R-l .0 ) /ECC ) A 262 

IF (G.LT.O.O) FT A=-FT A A 263 

AP=(AN-FTA)$CPR A 264 

FT A=FT A*OPR A 265 

ESTERR=TOL*T IME A 266 

WRITE ! 3 , 26 ) FT A , AP , TE » KERR , MQ , T I ME , ESTERR A 267 

26 FORMAT ( / 16H TRUE ANOMALY = F7 . 2 , 2X , 20HARG . OF PERIAPS1S = FIO.2,5 A 268 

1 X , 1 5HT0T AL EFFORT = F15.12//16H EXIT STATUS = I4,5X,15HN0. OF STE A 269 

2PS = 1 4 > 1 0 H TIME = ,FI5.10,14H EST. ERR = ,F8.6) A 270 

C IF THE INTEGRATION JUST PERFORMED WAS A SOLUTION REGENERATION A 271 

C THEN JUMP. A 272 

IF (KSW.EQ.l) GO TO 38 A 273 

C IF A BOUNDARY VALUE PROBLEM IS TO BE SOLVED USING THE FINAL A 274 

C TIME JUST COMPUTED, THEN JUMP. A 275 

IF (AF.NE.O.O) GO TO 27 A 276 

DD=NBS A 277 

PUNCH 36, ECCI,TAI,E,H,SMA, ECC, AN.FTA.AP, TIME, TE,DD A 278 

C A 279 

C AUTOMATIC INCREMENTATION OF INITIAL VARIABLES FOR ANOTHER A 280 

C INTEGRATION. A 281 

NCH-NCH- 1 A 282 

IF (NCH.LT.O) GO TO 1 A 283 

IF (NV.EQ.l) X I = X l +DV A 284 

IF (NV.EQ.2) XDOT I=XDOT I +DV A 285 

IF (NV.EQ.3) LR I = LR I +DV A 286 

IF (NV.EQ.4) F I = F I + DV A 287 

GO TO 4 A 288 

C A 289 

C WITH TIME ESTIMATE KNOWN, REGENERATE SOLUTION WITH FIXE.D A 290 

C STEP SIZE ANO STORE IN XOLD AS ESTIMATE OF SOLUTION TO A 291 

C THE BOUNDARY VALUE PROBLEM. A 292 

27 MQ=MQ*2 A 293 

IF (KERR.NE. I) GO TO 1 A 294 

IF (MQ.GT.1000) MC= 1000 A 295 

EST I M=T I ME A 296 

DEL T=EST I M/MQ A 297 

' DO 28 1 = 1 , 10 A 298 

28 Y ( 1 I =Y I N I T ( I I A 299 

XOLD 1 1 , 1 ) = Y I NIT ID A 300 

X0LD(2,i)=YINIT(2) A 301 

X0LD13,1)=YINIT13) A 302 

X0LD(4,1)=YINIT(5) A 303 

XOLD (5,1 )=YI NIT (6) A 304 

XOLD (6,1 ) = YI NIT t 7 I A 305 

TIME=0. 0 A 306 

M0DE=5 A 307 

DO 29 1=2, MQ A 308 

CALL RUNKUT ( T I ME , Y , P , NONL I N , NE , 0 .0, DELT , MODE ) A 309 

XOLD (1,1 )=U A 310 

XOLD ! 2,1 >=G A 311 

XOLO ( 3 , I ) =R A 312 

XOLD (4 , I ) =LU A 313 

XOLD ( 5 , I ) =LG A 314 

29 XOL D 1 6 , 1 ) = LR A 315 

C A 316 

C SET UP FOR QUASILINEARIZATION SOLUTION OF B. V. PROBLEM. A 317 

30 P ( I ) = MQ A 318 

P 1 2 ) =6 A 319 

P ( 3 I =3 A 320 

P ( 4 ) =0 A 321 

P 1 5 ) = 1 A 322 

P ( 6 ) = 1 A 323 

C A 324 

P( 19)=UE A 325 

P ( 20 ) =F A 326 

C SET UP FOR INTEGRATION OF MASS AS AN UNCOUPLED EQUATION. A 327 
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P(10)=l A 328 

P(ll)=l.O A 329 

C A 330 

C THE CONSTRAINTS ARE ENERGY, ANGULAR MOMENTUM, LAMBDA TRANS., A 331 

C AND HAMILTONIAN = 0 DUE TO FREE TIME. A 332 

CONST 11 ) =— 1 .0/(2. 0*AF) A 333 

CONST(2)=SQRT(AF*(1.0-ECCF**2) ) A 334 

CONST ( 3 ) =0. 0 A 335 

CONST (41=0.0 A 336 

CALL QUASI IP, NIT, PCT, CTO, ESTIM, CONST, DIFFI, BOUND) A 337 

IF (ESTIM. LT.O) 31,32 A 338 

31 GO TO 1 A 339 

C A 340 

C SAVE THE INITIAL VALUES FOUND BY QUASI FOR THE REGENERATION A 341 

C OF THE SOLUTION USING THE NONLINEAR EQUATIONS. A 342 

32 MSAVE=MQ A 343 

U=XOLD (1,1) A 344 

G=XOLD (2,1) A 345 

R=XOLD ( 3 , 1 ) A 346 

M= 1 • 0 A 347 

LU=XOLD (4,1) A 348 

LG=XOLD (5,1 ) A 349 

LR= XOLD (6,1) A 350 

LM= 0 . 0 A 351 

UEST=XOLD( 1 , MQ) A 352 

C • A 353 

C END POINT SUFFICIENCY CONDITION TEST. A 354 

IF (IF.NE.l ) GO TO 35 A 355 

P( 1 ) =MQ A 356 

NIB=NIT A 357 

CALL FOCAL ( P, NIB, PCT , CTO, ESTIM, CONST, DI FFI , FIXED, CONSTRT, ID.PARDE A 358 

1L) A 359 

IF (ID.EQ.2) WRITE (3,33) A 360 

IF (ID.NE.2) WRITE (3,34) ID A 361 

33 FORMAT 1//20X , *END POINT SUFFICIENCY CONDITION SATISFIED.*/) A 362 

34 FORMAT ( //20X , *END POINT SUFFICIENCY CONDITION NOT SATISFIED*, 13/ A 363 

2 20X ,45H********************************************* //) A 364 

C A 365 

C REGENERATE THE SOLUTION USING THE INITIAL VALUES FOUND BY A 366 

C QUASI AND THE NONLINEAR DIFFERENTIAL EQUATIONS. THE PURPOSE A 367 

C IS TO VERIFY THE NEWTON-R APHSON RESULTS AND GENERATE A A 368 

C COMPLETE LISTING OF THE TRAJECTORY. A 369 

35 X=ATAN(LG/(LU*U) ) A 370 

IF (LU.LT.O.O) X=X+Pt A 371 

P ( l ) =F A 372 

P ( 2 ) =UE A 373 

E=-.5 A 374 

AN=0.0 A 375 

TE=0. 0 A 376 

NB=l A 377 

NC=NCS A 378 

XDOT=-F*SIN(X)/(M*U!-LR*U*SINtX)*COS(G+X!/LG+U*SIN(X)*SlN(G+X)/R+C A 379 

1 OS ( G ) / ( U*R**2 ) A 380 

DD=Q. 0 A 381 

TA I =0. 0 A 382 

H=R*U*COS(G) A 383 

ECCI=SQRT(1.0-H*H! A 384 

CNO=NC A 385 

C A 386 

C PUNCH OUTPUT - INITIAL CONDITIONS A 387 

IF (ID.NE.2) GO TO 37 A 388 

PUNCH 36, F»UE|CN0»X , XDOT ,U*G,R,LU»LG,LR,LM A 389 

36 FORMAT (6E13.6) A 390 

C GO TO REGENERATE THE SOLUTION. A 391 

37 KSW=1 A 392 

GO TO 8 A 393 
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C A 394 

C PUNCH OUTPUT - FINAL CONDITIONS A 395 

38 IF ( I Do NE. 2 ) GO TO 39 A 396 

DD = NB A 397 

PUNCH 36, ECCI,TAI»E,H»SHA,ECC»AN,FTAjAP»TIME,TE,DD A 398 

A 399 

AUTOMATIC INCREMENTATION OF INITIAL VARIABLES TO SOLVE A A 400 

BOUNDARY VALUE PROBLEM NEAR THE ONE JUST GENERATED. A 401 

MQ=MS AVE A 402 

NCH=NCH— 1 A 403 

IF (NCH.LT.O) GO TO I A 404 

IF (NV.EQ.I) AF= AF+DV A 405 

IF (NV.EQ.2) ECCF=ECCF+DV ' A 406 

IF (NV.EQ.3) F= F -vD V A 407 

IF INV.E0.4) UE=UE+DV A 408 

GO TO 30 A 409 

END A 410- 

*$**$*>(<#***1}: B 1 

* COMPUTE PARAMETERS AND OUTPUT 4 B 2 

=4:2$:$ B 3 

B 4 
B 5 

SUBROUTINE OUTPUT B 6 

COMMON /MAIN/ E , H , SMA , ECC , CUT , S , APO , PER , X , W I ER , DPR , HAM B 7 

COMMON / I NT V AR/ T I ME , Y , P , NE , T OL , DX , MODE B 8 

REAL Y( 100) , P(20 ) B 9 

REAL M,LU,LG,LR,LM B 10 

EQUIVALENCE (Y(1),U), (Y(2),G), (Y(3),R), (Y(4),M), (YI5),LU), (Y( B 11 
16), LG), t Y ( 7 ) , LR ) , ( Y ( 8 ! , LM ) , (Y(9),AN), (Y(10),TE) B 12 

EQUIVALENCE (P(1),F), (P(2)*UE! B 13 

DATA PI /3. 141592653589793/ B 14 

E=.5*U*U-l.0/R B 16 

H=R*U*COS(G) B 16 

SMA=1.0/t2.0/R-U#U) B 17 

ECC=SQRT!1.0-H*H/SMA) B 18 

PER=SMA*( 1.0-ECC) B 19 

APO=SMA*( 1 . O+ECC ) B 20 

X=ATAN(LG/ (U*LU) ) B 21 

IF (LU.LT.O.O) X=X+PI B 22 

HAM=LU*(F*COS(X ) /M-SIN(G)/ ( R*R ) ) +LG<M F*SIN( X)/ ( M*U l-COSI G ) * ( 1.0/iR B 23 
1*R*U)-U/R) )+LR*U*SIN(G)-F*l 1.0+LM/UE) B 24 

W1ER=-F4LU/(M*C0S(X) ) B 25 

XD=X*DPR B 26 

GD=G*DPR B 27 

AND=AN*DPR B 28 

1 IF (XD.LT.180.) GO TO 2 B 29 

XD= XD-360. B 30 

GO TO l B 31 

2 IF (XD.GT.-180. ) GO TO 3 B 32 

XD=XD+360. B 33 

GO TO 2 B 34 

3 IF (ANO.LT.180. ) GO TO 4 B 35 

AND= AND-360 . B 36 

GO TO 3 B 37 

4 IF (AND. GT, -180. ) GO TO 5 B 38 

AND=AND*360. B 39 

GO TO 4 B 40 

5 WRITE (3,6) XO,E,H,U,GD,R,M,S,HAM,APO,AND,SMA,ECC,LU,LG,LR,LM,CUT, B 41 

1WIER.PER B 42 

6 FORMAT (1H , 8F 1 0 , 5 , F 10 . 7 , F 1 0 . 5/ 1H , 8 F 10 . 5 , F 1 0 . 7 , F 1 0. 5/ ) B 43 

RETURN B 44 

END B 45- 
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fc*## C 1 

❖ NORMAL-TANGENTIAL EQUATIONS OF MOTION* C 2 

Z = TIME, DX = STEP LENGTH, P = VECTOR OF PARAMETERS, C ^ 

Y = DEP. VAR., OER = VECTOR OF DERIVATIVES OF Y-S. C 5 

SUBROUTINE NONLIN (Z,0X,P) C 6 

REAL P120) ,Y ( 100! »DER( 100) C 7 

REAL M , LU, LG , LR, LM C 8 

COMMON /KUTTA/ DER, Y C 9 

EQUIVALENCE (Y(1),U), (Yi2),G), (YI3!,R), (Y(4),M), !Y(5>,LU), (Y( C 10 

16), LG), t Y ( 7 ) , LR ) , (Y(8!,LM), ( Y ( 9 ) , AN ) , (Y(10),TE) C 11 

C 12 

COMPUTE USEFUL COMBINATIONS. C 13 

F = P 1 1 ) C 14 

UE=P ( 2 ) C 15 

0PP=LG/U C 16 

HYP0T=1.0/SGRT10PP**2+LU**2) C 17 

XCOS=LL*HYPOT C 18 

XSI N=QPP*HYPOT C 19 

GCOS=COS ( G ) C 20 

GSIN-SIN(G) C 21 

FT=F*XC0S/M C 22 

FN=F*XSIN/(N*U) C 23 

RR= 1 . 0/ R C 24 

DT A=U*GCQS*RR C 25 

RR2=RR*RR C 26 

RUR2=RR2/U C 27 

C 28 

COMPUTE VELOCITY, FLIGHT PATH ANGLE, RADIUS AND MASS DERIVATIVES. C 29 

DER ( 1 )=FT-RR2*GS IN C 30 

DER(2)=FN-RUR2*GC0S+DTA C 31 

DER ! 3 ) =U*GS I N C 32 

DER (4 ) =~F/UE C 33 

C 34 

COMPUTE DERIVATIVES OF LAMBOA-S ASSOCIATED WITH ABOVE VARIABLES. C 35 

DERI5)=LG*(FN-RUR2*GCOS-DTA)/U-LR*GSIN C 36 

DER(6> = LU*RR2*GCOS-LG’MRUR2~U«‘RR)*GSIN-LR<‘U*GCOS C 37 

DER(7)=-2.0*LU*RR2*RR*GSIN-LG*12.0*RUR2*GC0S-DTA)*RR C 38 

DER(8)=(LU*FT+LG*FN)/M C 39 

C 40 

COMPUTE DERIVATIVES OF TRUE ANOMOLY, TOTAL EFFORT AND THRUST ANGLE C 41 

DER ( 9 ) =DTA C 42 

DER ( 10 ) =F C 43 

DER(U)=0.0 C 44 

DER(12)=0.0 C 45 

RETURN C 46 

END C 47- 

DRIVE A FUNCTION OF AN INTEGRAL TO ZERO D 2 

*$*$*##£#$##$#$###££$*#££$$£$*£$$$£#£ $£ 0 3 

D 4 

CONFUN = FUNCTION SUBR. WHICH COMPUTES THE VALUE OF THE FUNCTION D 5 

TO BE DRIVEN TO ZERO. D 6 

OLD = OLD VALUE OF THE FUNCTION BEFORE SIGN CHANGE. 0 7 

NEW = NEW VALUE OF THE FUNCTION AFTER SIGN CHANGE. D 8 

DXO = INTEGRATION STEP LENGTH THAT CAUSED FUNCTION VALUE TO D 9 

CHANGE FROM OLD TO NEW. D 10 

ZTOL = DESIRED MAXIMUM DEVIATION OF FUNCTION FROM ZERO. D 11 

DIFEQ = SUBROUTINE WHICH COMPUTES DERIVATIVES OF STATES DESCRIB- D 12 

ING THE SYSTEM. SEE INTEGRATION SUBROUTINE *RUNKUT*. D 13 

D 14 

SUBROUTINE CONVERG I CONFUN, OLD, NEW, DXO, ZTOL , D I FEQ ) D 15 

C 0 16 

COMMON /INTVAR/ T I ME , Y , P , NE , TOL , DUMMY , DUM D 17 

REAL Y(IOO) , P { 2 0 ) ,NEW D 18 

M0DE=5 D 19 

DX=DXO . D 20 

DO l 1=1,10 D 21 

DX=NEW*DX/iCLD-NEW) D 22 

OLD=NEW D 23 
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1 

2 


CALL RUNKUT ( T I HE , Y , P , D I FEQ , NE , T OL , DX , MODE ) 

NEW=CONFUN(T IKE, Y,P> 

IF (ABS(NEW) .LT. ZTOL ) 2,1 

CONTINUE 

RETURN 

OLD=0„0 

NEW=0.0 

RETURN 

ENO 

FUNCTION SWITCH (TIME.Y.P) 

REAL Y { 100) , PI20 ) 

SWI TCH=SQRT ( (Y(6 ) / Y ( 1 ) ) **2 + Y { 5 ) **2 )-Y I 4 ) * l 1. 0+Y ( 8 ) /P t 2 ) ) 

RETURN 

END 

FUNCTION CUTOFF (TIME,Y,P> 

REAL YC100)»P(20) 

CUTOFF = SQRT ( (Y(6)/Y( 1) ) **2*Y ( 5 ) **2 )- Y ( 4 ) 

RETURN 

END 

*** $**$>!:$$ *£$«:, >**>!.*<.$«:(:***** «<:*** 

QUASILINEARIZATION CONTROL SUBROUTINE * 

BY R. G. BRUSCH 11/25/68 

THIS SUBROUTINE CONTROLS THE GENERATION OF A SOLUTION TO A SET OF 
N NONLINEAR ORDINARY FIRST ORDER DIFFERENTIAL EQUATIONS WITH 
BOUNDARY CONDITIONS DESCRIBED AT THE FINAL AS WELL AS THE INITIAL 
POINT. THE SOLUTION TO THE NONLINEAR PROBLEM IS GENERATED AS THE 
LIMIT OF A SERIES OF SOLUTIONS TO A CONNECTED LINEAR PROBLEM. 

THE METHOD USED IS QUAS I L INEAR I Z AT I ON. 

REFERENCE, * A MODIFIED CUAS IL I NEARI Z AT ION METHOD FOR SOLVING 

TRAJECTORY OPTIMIZATION PROBLEMS*, JAY M. LEWALLEN, 
AIAA JOURNAL, VOL. 5, NO. 5, (MAY, 1967), PP 962-965. 

♦QUASILINEARIZATION AND NONLINEAR BOUNDARY VALUE 
PROBLEMS*, R. BELLMAN AND R. KAL ABA, 1965. 

♦SOLUTION CF VARIATIONAL PROBLEMS WITH BOUNDED CONTROL 
VARIABLES BY MEANS CF THE GENERALIZED NEWTON-R APH SON 
METHOD*, BY P. KENNETH AND G. E. TAYLOR, IN *RECENT 
ADVANCES IN OPTIMIZATION TECHNIQUES*, EDITED BY LAVI 
AND VOGL. 

DISCRIPTICN CF VARIABLES 
************************ 

PAR = A VECTOR CF PARAMETERS USED IN INTEGRATION. 1-fc RESERVED. 

THE FIRST b VALUES MUST CONTAIN THE FOLLOWING, INDEX = 

NS = NO. OF STEPS IN TIME CURRENTLY BEING USED. 1 

N = NUMBER OF FIRST ORDER D.E, (1ST M HAVE KNOWN I.C.) 2 

= NUMBER OF STATE VARIABLES 

M = NUMBER OF KNOWN STATE VARIABLE INITIAL CONDITIONS. 3 

IA = NUMBER OF *BANG-BANG* CONTROLS TIMES 3 4 

ITF = 0, IF FINAL TIME IS FIXED. = 1 IF FINAL TIME IS FREE. 5 

LD = MAX. NO. CF TIME STEPS FOR WHICH STATES CAN BE STORED. 6 


NIT = MAXIMUM NUMBER OF COMPLETE ITERATIONS ALLOWED FOR CONVERG. 
PC = THE MAX. PERCENT CHANGE IN INITIAL VALUES TO BE ALLOWED. 
THE PERCENTAGE IS BASED ON THE MAX. ABS. VALUE OF THE 
VARIABLE OVER ITS RANGE. HIGH PCT YIELDS RAPID CONVERG. 

LOW PCT. LESSENS PROBAB. OF UNCONTROL ABLE DIVERGENCE. 


D 24 
D 25 
D 26 
D 27 
D 28 
D 29 
D 30 
D 31 
D 32- 
E I 
E 2 
E 3 
E 4 
E 5- 
F 1 
F 2 
F 3 
F 4 
F 5- 
G 1 
G 2 
G 3 
G 4 
G 5 
G 6 
G 7 
G 8 
G 9 
G 10 
G 11 
G 12 
G 13 
G 14 
G 15 
G 16 
G 17 
G 18 
G 19 
G 20 
G 21 
G 22 
G 23 
G 24 
G 25 
G 26 
G 27 
G 28 
G 29 
G 30 
G 31 
G 32 
G 33 
G 34 
G 35 
G 36 
G 37 
G 38 
G 39 
G 40 
G 41 
G 42 
G 43 
G 44 
G 45 
G 46 
G 47 
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TOL = MAX VALUE OF THE SUM OF (MAX. DEVIATIONS OF THE STATES 


G 

68 

FROM THE STATES ON THE PREVIOUS ITERATION) 


G 

69 

USED TO TEST FOR CONVERGENCE COMPLETION. 


G 

50 

ESTIM = ESTIMATE OF FINAL TIME ( ASSUMMED FREE) 


G 

51 

CONST = A VECTOR CF CONSTANTS WHICH END POINT CONSTRAINTS MUST 

= . 

G 

52 

DIFFI = 


G 

53 

DIFFL = USER WRITTEN SUBROUTINES DEFINING LINEARIZED EQUATIONS 


G 

56 

BOUND = AND CONSTRAINT PROPERTIES. 


G 

55 

X = A VECTOR OF CURRENT STATE AND ALGEBRAIC VARIABLES 


G 

56 

(N STATE VARIABLES FOLLOWED BY IA ALGEBRAIC VARIABLES). 


G 

57 

AS SOLUTION IS DEVELOPED, SUCESSiVE VALUES ( WRT . TIME) 

ARE 

G 

58 

STORED IN XOLD MATRIX. 


G 

59 

XOLD = A MATRIX CF STATE AND ALGEBRAIC VARIABLES EVALUATED AT 

NS 

G 

60 

TIME INTERVALS WHICH PROVIDES THE SUBROUTINE WITH AN 


G 

61 

INITIAL SOLUTION WITH WHICH TO WORK. 


G 

62 



G 

63 

USE 


G 

66 



G 

65 

THE USER DESCRIBES HIS PARTICULAR PROBLEM BY SUPPLYING SUBR. DIFFI 

O 

66 

AND DIFFL, BY MODIFYING SUBR. BOUND, AND BY GENERATING A 


G 

67 

REASONABLE APPROXIMATION TO THE SOLUTION (XOLD) AND BY CALLING 


G 

68 

QUASI WITH THE PERTINENT SET OF PARAMETERS. 


G 

69 



G 

70 



G 

71 

SUBROUTINE QUASI ( PAR, HI T, PC, TOL , EST IM , CONST , DI FF I , BOUND ) 


G 

72 

«*##«*«*#**»$*** 


G 

73 

REAL XI 100) .CONST ( 10) ,UPDAT ( 10 ) , DER ( 100) .VARMAXtlO ) 


G 

76 

REAL XOLD (6, 1000) , PAR (20) ,ERR( 10), DUMMY! 100! .NEWER 


G 

75 

COMMON / IOLD/ XOLD 


G 

76 

COMMON /KUTTA/ DER, DUMMY 


G 

77 

EXTERNAL DLSUB 


G 

78 



G 

79 

SET UP CONSTANTS FOR THIS CALL. 


G 

80 

NS=PAR( 1 ) 


G 

81 

N=PAR(2 ) 


G 

82 

M=P AR ( 3 ) 


G 

83 

I A = PAR ( 6 ) 


G 

86 

I TF=P AR ( 5 ) 


G 

85 

NA=PAR 1 10) 


G 

86 

NVECT=N-M 


G 

87 

NVAR=(NVECT+1)*N+NA 


G 

88 

NUP=N-MvITF 


G 

89 

NST=N+ I A 


G 

90 

DO 1 1 = 1, NUP 


G 

91 

UPD AT ( I )=0.0 


G 

92 

UPDAT (N-M+l ) =0.0 


G 

93 

RH0=1 .0E+50 


G 

96 

RHOLD= 1 . 0E + 5 0 


G 

95 

CMAX= 1 • 0 


G 

96 

FRACT=0 . 0 


G 

97 

IDI V=0 


G 

98 



G 

99 

TEST FOR AN IMPROPER CALL. IF YES, STOP. 


G 

100 

IF (N.LT.2) STOP 701 


G 

101 

IF (M.LT.l) STOP 702 


G 

102 

IF (NS.LT.l) STOP 703 


G 

103 

IF (M.GE.N) STOP 705 


G 

106 

IF CNI T. LT. I } STOP 706 


G 

105 

IF (PC.LE.O.O) STOP 707 


G 

106 

IF (TOL. LT. 0.0) STOP 710 


G 

107 

IF (ESTIM. EC. 0.0) STOP 711 


G 

106 



G 

109 

SET LOOP TO ITERATE NO MORE THAN *NIT* COMPLETE TIMES. 


G 

110 

KLZ = 0 


G 

111 

DO 28 KRAP=KLZ , N I T 


G 

112 

IF (KRAP.EQ.O! GO TO 21 


G 

113 
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G 114 

INITIALIZE AL AND BL MATRICES FOR CONSTANT ELEMENTS. G 115 

CALL DIFFI (PAR) G 116 

G 117 

SET UP INITIAL STATE VECTOR FOR GENERATING SIMULTANEOUSLY G 118 

A PARTICULAR SOLUTION TO THE INHOMOGENEOUS LINEARIZED G 119 

EQUATIONS (USING XOLD(I.l) VECTOR AS THE INITIAL CONDITION) G 120 

AND N-M SETS OF SOLUTION TO THE HOMOGENEOUS LINEARIZED G 121 

EQUATIONS (THE INITIAL CONDITIONS FOR EACH SET OF HOMO. EQ. G 122 

BEING ALL ZEROS EXCEPT A ONE IN PLACE OF THE INITIAL STATE G 123 

WHICH IS FREE TO VARY.) • ' G 124 

G 125 

DO 2 1 = 1 ,NVAR G 126 

IF (I.LE.N) X ( I )=XOLD( 1,1) G 127 

IF (I.GT.N) X(I)=0.0 G 128 

CONTINUE G 129 

DO 3 I = 1 » NVECT G 130 

J=1*N+M+I G 131 

X(J)=1.0 G 132 

IF (NA.EQ.O) GO TO 5 G 133 

DO 4 1=1, NA G 134 

X1NVAR-NA+I )=PAR( 1 + 10) G 135 

G 136 

INITIALIZE ERROR ESTIMATING VECTOR AND VECTOR USED TO SAVE G 137 

THE MAXIMUM VALUE OF THE VARIABLES REQUIRING INI T. GUESSES. G 138 

DO 6 1=1, N G 139 

ERR ( I ) =0. 0 G 140 

DO 7 1 = 1 , NVECT G 141 

VARMAX ( I)=0.0 G 142 

TIME=0.0 G 143 

XNS=NS G 144 

DELT=EST 1 M/ XNS G 145 

G 146 

INTEGRATE ALL EQUATIONS FORWARD TO THE TIME ESTIMATE, ESTIM. G 147 

DO 10 K=2 , NS G 148 

PAR ( 1 ) = K-1 G 149 

G 150 

INTEGRATE EQUATIONS ONE STEP. G 151 

G 152 
G 153 

CALL RUNKUT ( T I ME , X , P AR , DLSUB, NV AR , 0 . 0, DELT , 5 ) G 154 

G 155 

SAVE MAXIMUM ERROR IN EACH DEPENDENT VARIABLE FROM PARTICULAR G 156 

SOLUTION AND STORE NEW PARTICULAR SOLUTION FOR USE ON NEXT ITER. G 157 

DO 8 1 = 1 ,N G 158 

NEWER=ABS(X( I )-XOLD( I , K ) ) G 159 

IF (NEWER. GT. ERR ( I ) ) ERR ( I ) =NEWER G 160 

XOLD ( I ,K)=X ( I ) G 161 

G 162 

SAVE MAX. ABSOLUTE VALUE OF VARIABLES INITIALLY GUESSED G 163 

DO 9 1=1, N G 164 

IF IABSIXII) >.GT.VARMAX( 1)1 VARMAX ( I ) =ABS ! X ( I ) ) G 165 

CONTINUE G 166 

0 CONTINUE G 167 

G 168 

COMPUTE ERROR METRIC G 169 

RH0=0.0 G 170 

DO 11 1=1, N G 171 

1 RHO=RHQ+ERK! 1 ) G 172 

G 173 

CHECK FOR HIGH DEVIATION OF INITIAL GUESS FROM EXACT NON- G 174 

LINEAR SOLUTION WITH SAME I.C. IF TRUE, DO NOT CHANGE G 175 

INITIAL VALUES. INTEGRATE ONCE TOWARDS NONLINEAR SOLUTION. G 176 

IF (KRAP.LE. 1. AND.RHO.GT. 10. *TOL ) GO TO 21 G 177 

G 178 

CHECK FOR A DIVERGING SOLUTION. ALLOW ONLY (5 CONSECUTIVE G 179 
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C DIVERGENCES. G 180 

IF t RHOLD. EQ.0.0 ) GO TO 13 G 181 

ER AT E= ( RHOLD— RHO ) /RHOL D G 182 

IF 1 ERATE.GT .0.0! IDIV = 0 G 183 

IF (ERATE.GT .0.0 ! GO TO 13 G 184 

I D I V= I DIV+1 G 185 

IF (IDIV.GE.15) WRITE (3,12) G 186 

12 FORMAT!//// 10X, * QUASILINEARIZATION CONTROL SUBR. UNCONTROLAB G 187 

2LE SYSTEM DIVERGENCE.*////! G 188 

IF (IDIV.GE.15) GO TO 30 G 189 

C G 190 

C COMPUTE CHANGES IN INITIAL CONDITIONS WHICH WOULD SATISFY G 191 

C B. C. IF THE SAME XOLD ESTIMATE WERE USED AGAIN. G 192 

13 CALL BOUND ( X , DER , CONST , PAR , UPDAT ) G 193 

C G 1 94 

C HOWEVER, ONLY USE A FRACTION OF THESE INITIAL VALUE UPDATES. G 195 

C THIS WILL PRESERVE CONVERGENCE EVEN WITH BAD INIT. GUESSES. G 196 

C G 197 

C IF SOLUTION IS CONVERGING MAKE ADJUSTMENTS TO INITIAL STATES G 198 

C AND TIME ESTIMATE PROPORTIONAL TO FRACT. G 199 

C FRACT IS CHOOSEN SO THAT NO CORRECTION WILL EXCEED A PER G 200 

C CENT OF THE MAX. ABSOLUTE VALUE OF THAT VARIABLE OVER RANGE. G 201 

CMAX=1 . 0E-50 G 202 

DO 14 1=1, NVECT G 203 

J=I*M G 204 

IF (ABSIUPDATt I )/VARMAX(J) ) .GT.CMAX) CMAX= AB S ( UPDAT ( I ) /VARMAX ( J ) ) G205 

14 CONTINUE G 206 

IF IITF.EQ.l ) 15,16 G 207 

15 CMAX=AMAX1(C MAX, ABSIUPDATt NVECT+1)/ESTIM)J G208 

C G 209 

C IF PERCENT ERROR IN LAST ITERATION .GT. PCT, REDUCE FRACT. G 210 

16 DO 17 1=1, N G 211 

IF (ABS(ERR( I )/VARMAX( I ) I .GT.CMAX) CKAX= ABS t ERR ( I ) /VARMAX ( I ) ) G212 

17 CONTINUE G 213 

PCT=PC*0.01 G 214 

IF (CMAX.GT. PCT ) FRACT = PCT/CM,AX G 215 

IF (CMAX.LT.PCT) FRACT=1.0 G 216 

IF (KRAP.GT. NIT-10. AND. FRACT. EQ. 1.0) FRACT=0.8 G 217 

C G 218 

C UPDATE INITIAL VARIABLES. G 219 

RHOLD= RHO G 220 

IF (ITF.EQ. 1 ) 18,19 G 221 

18 ESTIM=ESTIM+UPOAT(NUP)*FRACT G 222 

19 DO 20 1 = 1, NVECT G 223 

K=M+I G 224 

20 XOLDIK, l)=XOLD(K , 1 ) +UPDATI I )*FRACT G 225 

C G 226 

C TEST OUTPUT. G 227 

21 ND = NS/9 + 1 G 228 

K=N + 1 A G 229 

IF (KRAP.EQ.O) WRITE (3,22) G 230 

22 FORMAT (1HU G 231 

WRITE (3,23) KRAP, RHO, FRACT, EST I M , UP D AT ( N-M* 1 ) G 232 

23 FORMAT!/* ITERATION = *, 15, * RHO = *, E12.3, * FRACTION = *, G 233 

2 F 6 . 4 , * TIME EST. = *, E16.8,* TIME UPDATE = *,E12.4) G 234 

DO 27 1=1, K G 235 

IF ( I.LE.M.OR. I .GT.N) 26,24 G 236 

24 WRITE (3,25) ( XOLD ( I , J ) , J= 1 »NS ,ND ) , XOLD ( I , N S ) , UPDAT ! I -M) G237 

25 FORMAT (4X , 1 E16. 8 , 8E10. 2 , 1E16. 8, IE 10 . 2 ) G238 

GO TO 27 G 239 

26 WRITE (3,25) ( XOLD( I , J ) , J= 1 , NS , ND ) , XOLDt I , NS ) G240 

27 CONTINUE G 241 

C G 242 

C CHECK FOR METRIC WITHIN TOLERANCE. IF NOT REPEAT. G 243 

IF (RHO. LT.TCL. AND. CMAX.LT.TOL. AND. KRAP.GT. 2) GO TO 31 G 244 

28 CONTINUE G 245 

C G 246 

WRITE (3,29) G 247 

29 FORMAT (////10X, *CU AS I L I NE AR I Z AT I ON CONTROL SUBR. THE SYSTEM W G 248 

20ULD NOT CONVERGE IN THE ALLOTTED NUMBER OF ITERATIONS.*////) G 249 

30 EST I M=-l . 0 G 250 

31 PAR ! 1 ) =NS G 251 

RETURN G 252 

END G 253- 



non won ro *- o o o o ooooooooooonooooooo 


205 


* FOURTH DEGREE RUNGE KUTTA INTEGRATION * H ‘ 2 

H 4 

BY R. Go BRUSCH 2/13/68 H 5 

MINIMUM TRUNCATION ERROR AS PER, H 6 

Fo CESCHINO, NUMERICAL SOLUTION OF INITIAL VALUED PROBLEMS, H 7 

PAGES A A , 45 , AND 67. H 8 

H 9 

MODE = 5 FOR A FIXED STEP LENGTH = DX H 10 

= 0 TO START INTEGRATION WITH STEP LENGTH SELECTED EVERY H 11 

4TH INTEGRATION SO ERR/UNIT TIME IS LESS THAN THE TOL . H 12 

NOTE THAT THIS DOES NOT NECESARILY GUARANTEE MINIMUM TOTAL ERROR H 13 

FOR THE NUMBER OF STEPS USED, SINCE ERROR PROPAGATION IS NOT H 14 

TAKEN INTO C ONS I DER AT I CN . H 15 

MODE SHOULD NOT BE ADJUSTED BY THE MAIN LINE PROGRAM IN H 16 

THE VARIABLE STEP LENGTH MODE. H 17 

H 18 
H 19 

SUBROUTINE RUNKUT ( X , Y , P , D I F , NDE , TOL , DX , MODE ) H 20 

REAL X,Y( 100) ,P(20) ,Y1( 100) , DERI 100) , XK (4, 100) , TOL, A(4,4) ,YOLD( 3, 1 H 21 

100) ,DER0LD(3,100) ,XS(4) H 22 

COMMON /KUTTA/ DER , Y 1 H 23 

DATA A( 1 , 1 ), A(2, 1 ) ,A(2, 2 1/0.40,-. 14999999999999,0. 750/, A (3, 1 ) ,A(3, H 24 

12) ,A(3, 3)/. 431818181818181, -.340909090909090,0. 909090909090910/, A( H 25 

24,1) ,A( 4,2) ,A(4,3), A (4,4 I/O. 152777777777777, .34722222 2222221,-34 72 H 26 

322222222221,-152777777777777/ H 27 

H 28 

CHECK MODE. IF INITIALLY ZERO, TAKE SMALL STEPS 4 TIMES TO SET H 29 

STEP LENGTH COMPATIBLE WITH SINGLE STEP ERROR TOLERANCE. H 30 

IF INITIALLY 5, A FIXED STEP LENGTH = DX IS USED. H 31 

MQ=0 H 32 

IF ( MODE . GT . 0 ) GO TO 2 H 33 

DX= . 000 1 H 34 

H = • 000 1 H 35 

HN=1.0E+50 H 36 

MODE=MODE-l H 37 

GO TO 3 H 38 

IF ( MODE. EQ. 5 ) H = DX H 39 

IF ( MODE .GT . 5 ) STOP H 40 

H 41 

SET THE INDEPENDENT VARIABLES H 42 

XS(1)=X H 43 

XS(2!=X+H*0.4 H 44 

XSI3)=X4-H4.6 H 45 

XS ( 4 ) =X + H H 46 

L= I ABS i MODE ) H 47 

H 48 

INTEGRATE - DIF IS A SUBROUTINE WHICH COMPUTES THE D/DX-S USING H 49 

Yl-S, X, H, AND PARAMETERS P H 50 

DO 4 J = I , N D E H 51 

4 Y1(J)=Y(J) H 52 

DO 8 1=1,4 H 53 

CALL DIF (XSm.H.P) H 54 

DO 8 J= L , NDE H 55 

XKt I , J)=DER( J) H 56 

DER ! J ) = 0. 0 H 57 

DO 5 K=l,l H 58 



noooooooooooooooonooonooonoooo os ■**! onoo oo 


206 


5 DER(J)=DER(J)+A( I,K)*XK(K,J) H 59 

IF (I.LT.4) Y 1 ( J)=Y( J)+H*DER( J ) H 60 

IF (I.LT.4) GO TO 8 H 61 

Y( J)=Y( J)-t-H*DER(J> H 62 

H 63 

SAVE Y-S AND THEIR DERIVATIVES FOR ERROR EVALUATION AND H ADJUST. H 64 

IF (L.EQ.5) GO TO 8 H 65 

IF (L.EQ.4) GO TO 6 H 66 

YOLDIL, J)=Y( J) H 67 

DEROLDIL, J)=DER( J) H 68 

GO TO 8 H 69 

H 70 

COMPUTE NEW H TO MAKE MAXIMUM ERROR LESS THAN THE TOLERANCE. H 71 

FOR 4TH ORDER RUNGE KUTTA, UNIT ERROR IS PROPORTIONAL TO H**4. H 72 

STEPER=+( . 18333333333333*(Y (J )-YOLD( 1, J ) ) + . 4 5*( YOLD ( 3 , J ) - YOLD 1 2 » J ) H 73 

1 1 )-( .05*(D£R ( J) +DEROLDI 1 , J) ) + . 45* ( DEROLD ( 3 , J ) +DER0LDI 2, Ji ) ) *H H 74 

IF (STEPER.EQ.0.0) GO TO 7 H 75 

UNI TER=STEPER/H H 76 

HNE'W=H# ( SORT ( ABS ( TOL/UNI TER ) ) ) H 77 

IF (ABS(HNEW).LT.HN) HN= ABS t HNEW ) H 78 

IF (ABS (HNEW) .LT.HN) MC=J H 79 

IF (ABS(HNEW). LT.HN) UN=UN I TER H 80 

IF (J.LT.NDE) GO TO 8 H 81 

X= X + H-HN H 82 

H=HN H 83 

DX=H H 84 

HN=1.0E+50 H 85 

M0DE=0 H 86 

CONTINUE H 87 

X = X+H H 88 

IF ( MODE . LT . 0 ) GO TO 1 H 89 

IF ( MODE . LT . 4 ) MODE = MODf. + l H 90 

RETURN H 91 

END H 92- 

****************************************************************** I 1 
COMPUTE DERIVATIVES FOR PARTICULAR SOLUTION AND N-M HOMO. SOLUTION I 2 

****************************************** **** ************* ******* J 3 

I 4 
1 5 

THE *N* LINEARIZED DIFFERENTIAL EQUATIONS ARE OF THE FORM, 1 6 

D(X(J))/DT = AL(J,!)*X(I) + BL ( I ) (1) I 7 

WHERE REPEATED SUBSCRIPTS IMPLIES SUMMATION. I 8 

X ( I ) = THE NEW ESTIMATE OF THE DEPENDENT VARIABLES. 1 9 

ALIJ, I) = A MATRIX OF COEFFICIENTS WHICH ARE EVALUATED USING ONLY l 10 

THE OLD STORED VALUES OF THE DEPENDENT VARIABLES. I 11 

THEY REPRESENT THE THE FIRST PARTIAL DERIVATIVES OF THE I 12 

RIGHT HAND SIDES OF THE NON-LINEAR D. E. WRT . THE STATE 1 13 

BL(J) = A VECTOR OF CONSTANT COEF WHICH ARE EVALUATED USING ONLY I 14 

THE OLD STORED VALUES OF THE DEPENDENT VARIABLES. I 15 

I 16 
I 17 

THE MATRIX AL AND VECTOR BL ARE FOUND FROM THE NON-LINE I 18 

•EQUATIONS BY A TAYLOR SERIES EXPANSION OF THE RIGHT HAND SIDE OF T I 19 
1-ST ORDER NONLINEAR DIFFERENTIAL EQUATIONS ABOUT THE OLD STORED I 20 
SOLUTION AND RETAINS ONLY LINEAR TERMS OF THAT EXPANSION. I 21 

SINCE THE RESULTING EQUATIONS ARE LINEAR IN THE DEPENDENT VARIABLE I 22 
THE BOUNDARY CONDITIONS CAN BE SATISFIED IDENTICALLY AT EACH STEP I 23 
BY THE PRINCIPLE OF SUPERPOSITION. I 24 

THIS SUBROUTINE GETS THE VALUES FOR AL AND BL ONCE AND I 25 

USES THEM TO GENERATE THE DERIVATIVES OF THE INHOMOGENEOUS SOLUTIO 1 26 

VECTOR AS WELL AS N-M HOMOGENEOUS SOLUTION VECTORS WHICH WILL BE I 27 
NEED FOR THE SUPERPOSITION. THUS THE A AND B MATRICES ONLY I 28 

HAVE TO BE EVALUATED ONCE INSTEAD OF N-M+l TIMES IF ALL OF THE I 29 

SOLUTIONS WERE NOT BEING GENERATED SIMULTANEOUSLY. I 30 

SUBROUTINE DLSUB ( T , DELT , PAR ) I 31 

C **************** 132 
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REAL T,X( 100 ) ,DER( 100) ,PAR( 20! , AL (6, 6) , BU 6) I 33 

COMMON /I NT EG/ AL.BL l 34 

COMMON /KUTTA/ DER.X I 35 

IS=PAR ( 1 ) I 36 

N=PAR ( 2 ) I 37 

M=P AR ( 3 ) l 38 

NVECT=N-M I 39 

I 40 

GET CURRENT VALUE FOR AL, EL MATRICES. I 4L 

CALL DIFFL I IS, PAR) I 42 

I 43 

GET DERIVATIVES OF INHOMOGENEOUS EQUATIONS I 44 

DO 2 1=1, N I 45 

SUM=0.0 I 46 

DO 1 J=1,N I 47 

SUM=SUM«-AL(I,J)*X(J) I 48 

DER ( I ) = SUM + BL ( I ) I 49 

I 50 

GET DERIVATIVES OF <‘NVECT* SETS OF HOMOGENEOUS EQUATIONS. I 51 

DO 5 K=1,NVECT I 52 

KV=K*N I 53 

DO 4 1=1, N I 54 

SUM=0. 0 I 55 

DO 3 J=L,N I 56 

SUM=SUM+AU I , J)*X( J+KV) I 57 

DER(I+KV)=SUM I 58 

CONTINUE I 59 

I 60 

USER INSERTS DERIVATIVES OF *PAR(10>* NO. OF AUXILLARY I 61 

VARIABLES WHICH ARE UNCOUPLED FROM THOSE INVOLVED IN THE I 62 

QUASILINEARIZATION. THE FIRST IS CALLED X II N-Mf.1 )*N+ 1 ) I 63 

DER ( 2 5 ) =-P A R (201/PARI 19) I 64 

RETURN I 65 

END I 66- 

J l 

THE LINEARIZED DIFFERENTIAL EQUATIONS * J 2 

4## 4###$ 4' J 3 

J 4 

USER WRITTEN ROUTINE TO EVALUATE AL AND BL MATRIX USING THE OLD J 5 

STORED VALUES. J 6 

J 7 

PAR = A VECTOR OF CONSTANT PARAMETERS NEEDED FOR DERIVITIVE. J 8 

XOLDII.J) = THE VALUE OF THE I-TH STATE AFTER J INTEGRATION J 9 

STEPS ( FROM LAST ESTIMATE OF THE SOLUTION) J 10 

J 11 
J 12 

SUBROUTINE DIFFL (IS, PAR) J 13 

COMMON /INTEG/ AL , BL J 14 

COMMON / 1 OLD / XOLD J 15 

COMMON /KUTTA/ DER, Y J 16 

REAL ALI6.6 ) ,BL (6) , PARI 20) ,XOLD( 6, 1000) , DERI 100) , Y( 100) J 17 

REAL M, LU , LR , LG J 18 

M=Y(25) J 19 

UE= PAR (19) J 20 

F=P AR ! 20 ) J 21 

U=XOLDtl,IS) J 22 

G=X0LD(2,IS) J 23 

R=X0LD(3,IS) J 24 

LU= XOLD (4,15) J 25 

LG= XOLO ! 5 , 1 S ) J 26 

LR=XO' 0(6, IS) J 27 

COSG=COS ( G ) J 28 

J 29 

COMPUTE FREQUENTLY APPEARING FACTORS. J 30 

SI NG=S I N t G ) J 31 

FOM=F/M J 32 
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D=(LG*LG)/(U<‘U)-t-LU*lU J 33 

SQRTD=SQRT(D) J 34 

D32=SQRTD*D J 35 

J 36 

UDOT EQUATION. J 37 

AL(1»1)=F0N*LU*LG**2/(D32*U**3> J 38 

AL(1,2)=-C0SG/(R**2) J 39 

AL(lj3!=2.0*SING/R**3 J AO 

AL(1,4)=F0M*!LG/U>**2/D32 J 41 

ALU > 5)=-F0M*LU*LG/(D32*U**2) J 42 

J A3 

G DOT EQUATION. J A4 

ALU, l )=-F0R*LG*(LG>i'*2 + 2.0*(LU*U)**2 ) / ! D32*U**5 > + ( COSG/ ( (R*U)**2)4- J A5 

1C0SG/R) J A6 

AL(2,2)=SING/(U*R**2>-U*SING/R J A7 

AL(2,3!=(2.0/(U*R**3)-U/R**2!*COSG J A8 

AL [2,4>=-F0f'*LG*LU/ (D32*U**2) J A9 

AL(2,5)=F0M*(LU/U)**2/D32 J 50 

J 51 

ROOT EQUATION. J 52 

AL(3jl ) = S I NG J 53 

AL ( 3 , 2 ! =U*C0SG J 5A 

J 55 

LU OOT EQUATION. J 56 

AL t A, 1 !=-F0M*(LG**2 I * ( 2 . 0 « ( LG/U ) **2 + 3 . 0*LU** 2 ) / ( D32*U**4 ) + 2 . 0*LG *C J 57 

10SG/(U*(U*R) **2) J 58 

AL(4,2)=LG*t 1./ IR*U)**2+1.0/R)*S ING-LR*C0SG J 59 

AL(4,3!=LG*(2./tR*tR*U)**2)+l.0/R**2)*C0SG J 60 

AL(4,4)=-F0f'*LG**2*LU/ (D32*U**3) J 61 

AL ( A, 5 ) =F0M*LG* ( ( LG/U) **2+2.0*LU**2 ) t ( D32*U**3 )- I 1 . 0/ ( R*U ) **2*1 . 0/ J 62 

1R)*C0SG J 63 

AL(4,6)=-SING J 6A 

J 65 

LG DOT EQUATION. J 66 

AU5,l> = LG*U.0/(lR*U>**2) + l./R>*SING-LR*C0SG J 67 

AL(5»2!=(LR*U-LU/R**2>*SING-LG*( 1 .0/ ( R*R*U 1-U/R )*COSG J 68 

AL(5,3)=-2.0*LU*C0SG/R**3tLG*(2.0/(U*R**3>-U/R**2)*SING J 69 

AL(5,4)=C0SC/R**2 J 70 

AU5,5)=-(1.0/IU*R**2)-U/R)*SING J 71 

AL(5)6!=-U*C0SG J 72 

J 73 

LR DOT EQUATION. J 7A 

AL(6,1)=LC*<2.0/(R*(R*U)**2)+1.0/R**2)*COSG J 75 

AL ( 6 , 2 ) = AL ( 5 , 3 ) J 76 

ALI6 I 3)=6.0*LU*SING/R**4+LG*(6.0/(U*R**4)-2.*U/R**3)*C0SG J 77 

AL16,4)=-2.0*SING/R**3 J 78 

ALl6,5!=-(2.0/(U*R**3)-U/R**2}*C0SG J 79 

J 80 

NON-UNEAR TERNS. J 81 

BL (1 I=F0M*LU**3/032-3.0*SING/R**2+G*C0SG/R**2 J 82 

8L(2!=F0M*LG*(2.*{LG/U)**2+3.0*LU**2)/(D32*U**2 1+1-4. 0/tU*R**2)+U/ J 83 

1R)*C0SG-G*AL(2,2) J 8A 

BL(3)=-U*G*CCSG . J 85 

BL (4)=FOM*LG**2* (2.0*! LG/U > **2+3.0*LU**2 ) / ( D 32*U**3 ) +LG* !-A. 0/(R*U J 86 

1)**2-1.0/R)*C0SG-G*AL(4,21 J 87 

BL(5)=-3.0*LG*SING/ I U*R**2 ) + ( 2 .0*LU /R**2+LR*U ) *COSG-G*AL ( 5 » 2 ) J 88 

Bl{6)=-LG*<8./(U*R**3I-U/R**2)*COSG-6.*LU*SING/R**3-G*AU6»2) J 89 

RETURN J 90 

END J 91- 
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INITIALIZE AL AND SI- MATRICES FOR CONST*. ELEMENTS K 2 

&$$$$ £##$$$$ $$$$$$$£$$$$$$$$$$$$$£$$ $ )fit # $= $ # -If. >Jt # U £ $ ❖ *C=# =( s ❖ fc fc $ K 3 

K A 
K 5 

SUBROUTINE OIFFI (PAR) K 6 

COMMON / 1 OLD / XOLD K 7 

COMMON /INTEG/ AL.BL K 8 

REAL \L(6,6> ,BL I 6 ) , PAR ( 20 > , XOLD! 6 , 1000) K 9 

N=PAR ( 2 ) K 10 

M=PAR ( 3 ) K 11 

K 12 

DO 1 1 = 1, N K 13 

BL ( I > =0.0 K 14 

DO 1 J=l,N K 15 

AL! I , J ) =0.0 K 16 

K 17 

RETURN K 18 

END K 19- 

$$£##9 L 1 

END POINT SUFFICIENCY CONDITION TEST * L 2 

if fc* L 3 

L A 
L 5 

SUBROUTINE FOCAL ( PAR, NI T , PCT, TOL , EST IM , CONST, DI FF I , F IXED.CONSTRT , L 6 

1 1 D , D) L 7 

REAL X0LDI6, 1000 ), PARI 20 ), CONST! 10) , SAV( 10) L 8 

REAL PARTWOI 10, 10 i ,PHI ( 10, 10) , A! 10, 10) ,XS( 10) L 9 

INTEGER I CON ST (10) L 10 

COMMON /10LD/ XOLD L 11 

NS = PAR ( 1 ) L 12 

N= P AR ( 2 ) L 13 

M=P AR ( 3 ) L IA 

I TF =.P AR ( 5 ) L 15 

NTC = N-MMTF L 16 

L 17 

GET NO. OF VARIABLES INVOLVED IN CONSTRAINTS, NO. OF L 18 

INDEPENDENT VARIABLES INVOLVED IN CONSTRAINTS, LOCATIONS OF L 19 

EACH, 2ND PARTIAL DERIVATIVES OF G, AND PHI CONSTRAINT MATRIX L 20 

CALL CONSTRT ( L , P AR , NV , N I , I CONST , P ART WO, PH I ) L 21 

L 22 

GENERATE PARTIAL DERIVATIVES OF LAMBDA FINALS WRT. DEPENDENT L 23 

AND INDEPENDENT VARIABLES INVOLVED IN CONSTRAINTS. L 24 

WRITE (3,1) ( (PARTWOI J, I ) , 1 = 1, NV) , J= l,NV) L 25 

FORMAT (3E16.8) L 26 

WRITE (3,2) ( (PHI ( I , J) , J=1,NI ) , I =1,NV) L 27 

FORMAT ( / ( E 1 6. 8 ) ) L 28 

DO 3 I = 1 , N L 29 

SAVin = XOLb( I, NS) L 30 

TI MS AV= E ST I M L 31 

DO 4 1=1, N L 32 

XS 1 I ) = XOLD (1,1) L 33 

IF (NV.GT.NTC) STOP 40 L 34 

DO 10 1=1, NV L 35 

DO 7 J =1 , NV L 36 

CONST ( J ) =S AV ( J ) L 37 

IF (I.EQ.J) 5,7 L 38 

CONST! J)=SAV(J)*( i.O+D) L 39 

DENOM=SAV( J) *D L 40 

L 41 

TEST FGR ENDPOINTS NEAR ZERO. L 42 

IF ( ABSiCONST 1 J ) ) .LT. ,001 ) 6,7 L 43 

6 CONST! J)=SAVIJ)+. 01 L 44 

DENOM=.Ol L 45 

7 CONTINUE L 46 

CALL QUASI (PAR, NIT, PCT, TOL , EST I M, CONST , D IFF I , FIXED) L 47 

IF (ESTIM.LT.O.O) I D=0 L 48 

IF (ESTIM.LT.O.O! RETURN L 49 

NH=N/2 L 50 

K=0 L 51 

DO 9 L= 1 , NH L 52 

IF ! I CONST I L 1 . NE . 0 ) 8,9 L 53 

K=K+1 L 54 
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PART WO (K,I)=PARTKC(K,I ) -M XOLD ( L + NH, N S ) -S A V ( L + NH ) 1/DEN0M L 55 

9 CONTINUE L 56 

IF IK.GT.NV) STOP 62 L 57 

EST I M=T I MSAV L 58 

DO 10 K=l,N L 59 

10 XOLD ( K , 1 ) =XS ( K ) L 60 

WRITE (3,1) ( (PARTWOt J, I ), 1=1, NV ) , J=l,NV) L 61 

L 62 

POST- AND PRE-MULTIPLY BY PHI L 63 

DO 12 1=1, NV L 64 

DO 12 J= 1 > N I L 65 

SUM=0.0 L 66 

DO 11 K= 1 , N V L 67 

1 SUM=SUM+PARTWO( I ,K)#PHI (K, J) L 68 

2 A ( I » J) = SUM L 69 

DO 14 1=1, NI L 70 

DO 14 J= 1 ,N I L 71 

SUM=0.0 L 72 

DO 13 K= 1 , NV L 73 

3 SUM=SUM+PH1 ( K, I ) #A(K, J ) L 74 

4 PART WO ( I , J i = SUM L 75 

WRITE (3,2) UPARTWOI J, I ) , 1=1, NI ) , J=1,NI ) L 76 

L 77 

CHECK DEFINITENESS. L 78 

CALL DEFINIT ( P ARTWO , N I , I D ) L 79 

RETURN L 80 

END L 81- 

6#666*6#66#666*66#666#666666*666#66#6666 M 1 

COMPUTATION OF REAL VARIABLE CONSTRAINTS * M2 

*##66 #666# 6 #66# *#666 6 6 66666 6 666666666666 M 3 

KS =1 FOR A NORMAL CALL FROM FOCAL M 4 

PAR = A VECTOR OF PARAMETERS AS DESCRIBED IN QUASI. M 5 

NVAR = THE NUMBER OF VARIABLES INVOLVED IN THE CONSTRAINTS. M 6 

NI = THE NUMBER OF INDEPENDENT VARIABLES INVOLVED IN CONSTRAINTS M 7 

NC = NO. OF CONSTRAINED VARIABLES = NO. OF CONSTRAINTS. M 8 

1 CONST = 1 IF X(J) IS A DEPENDENT VARIABLE IN A CONSTRAINT. M 9 

= 2 IF X(J) IS AN UNCONSTRAINED VARIABLE IN CONSTRT. M 10 

= 0 OTHERWISE M 1 1 

PARTWOI I , J ) = VALUE OF SECOND PARTIAL DERIVATIVE OF G WITH RESPECT M 12 

TO THE I-TH AND J-TH VARIABLE INVOLVED IN THE M 13 

CONSTRAINTS. THE MATRIX IS ORDERED SO THAT ALL M 14 

DERIVATIVES WRT THE DEPENDENT VARIABLES APPEAR IN M 15 

THE UPPER-LEFT. (WITH I AND J .LT. NC) M 16 

PHI ( I , J ) = A MATRIX RELATING CONSTRAINED VARIABLES TO THE M 17 

UNCONSTRAINED VARIABLES. CON ( I ) = PH It I , J ) *UNC ON ( J ) M 18 

M 19 
M 20 

SUBROUTINE CONSTRT ( KS , P AR , NVAR, N I , ICONST , PART WO, PHI ) M 21 

*6#666##6###6*66*# M 22 

REAL PAR120) , PARTWOI 10, 10) , PHI ( 10, 10 ) ,X0LD(6, 1000) ,PARCD ( 10, 10) M 23 

REAL PARCI 1 10, 10) .PARCDI ( 10, 10), X( 10) M 24 

INTEGER I CONST ( 10 ) M 25 

COMMON / 1 OLD/ XOLD M 26 

REAL MU 1 , MU 2 M 27 

REAL LG , LR M 28 

NS = P AR 1 1 ) M 29 

N=P AR ( 2 ) M 30 

USER LOADS ICONST, PARTWO, PARCD, AND PARCI HERE. M 31 

#6*6*66 66 4 #6 666 6 6**6 6 6 #*#6**66 66 666# *6 #6*666 66*6# M 32 

M 33 

CHOOSE U AND G TO BE INDEPENDENT M 34 

U=XOLD( 1 ,NS ) M 35 

G=XOLD ( 2 , NS ) M 36 

R=XOLD ( 3 , NS ) M 37 

LG=XOLD ! 5 , NS ) M 38 

LR=XOLD ( 6 ,NS ) M 39 
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MU1=-R*R*(LG/(R*TAN(G)!+LR! M 40 

MU2=LG/(R*U*SIN(G) ! M 41 

PART WO I 1 j 1 ) = MUl H 42 

PAR TWO { 2 » 2 ) = —MU 2 *l/*R*COS t G } M 43 

PARTW0(3,3)=-2.0*MU1/R#*3 H 44 

PARTW0(1,2)=PARTW0(2,1)=-MU2*R*SIN(G) M 45 

PART WO 1 1 , 3 ) = PAR TWO I 3, 1 ) = MU2*C0S( G ) M 46 

PART WO (2,3)=PARTWO(3»2 ) =— MU2 *U*S I N ( G ) M 47 

ICONST ( 1 )=1 M 48 

1 CONST ( 2 ) = 1 M 49 

I CONST ( 3 ) =2 M 50 

ICONST ( 4 ) =0 H 51 

I CONST ( 5 ) = 0 M 52 

ICONST (61 = 0 H 53 

H 54 

COMPUTE PARTIAL DERIVATIVES OF CONSTRAINTS WRT. DEPENDENT VAR M 55 

PARCD ( 1,1) = U M 56 

PARCD (1,21=0.0 M 57 

PARCD(2,1)=R*C0S(G) M 58 

PARCD(2,2)=-R*U*SINIG> M 59 

M 60 

COMPUTE PART 1 AL S OF CONSTRAINTS WRT. INDEPENDENT VARIABLES. M 61 

PARCI (1 , 1 )=1 .0/R**2 M 62 

PARCI 12,1)=U*CQSIG) M 63 

fc £#£4 # ❖t' $$$$*$:(<$❖ 4 4$ M 64 

M 65 

SET NC, NVAR, AND NI. M 66 

NC=0 M 67 

NVAR=0 M 68 

NI =0 M 69 

DO 1 1=1, N M 70 

IF ( I CONST (I ) . NE .0 ) NV AR=NV AR+ 1 M 71 

IF I ICONST II ) • EQ. I ) NC = NC + 1 M 72 

IF t ICONST ( I ) .EC. 2) NI=N1+1 M 73 

CONTINUE K 74 

M 75 

COMPUTE PHI = -PARCD! INVERSE )*PARC I M 76 

KC = NC M 77 

CALL SIMEQ ( PARCD , X , KC , PARCD I , X ) M 78 

IF ( KC . NE. NC ) STOP 31 M 79 

DO 3 1=1, NC M 80 

DO 3 J = 1 , N I M 81 

SUM=0 . 0 M 82 

DO 2 K=1,NC M 83 

SUM=SUM-PARCDI ( I ,K)*PARCI (K, J) M 84 

PHI(I,J)=SUM M 85 

M 86 

FILL OUT LAST ROWS OF PHIfl.J) WITH IDENTITY MATRIX. M 87 

DO 4 1 = 1 ,NI M 88 

DO 4 J= 1 , N I M 89 

PHI C NC ■*- 1 , J i =0.0 M 90 

IF II.EQ.J) PHI (NC+I ,J 1=1.0 M 91 

4 CONTINUE M 92 

RETURN M 93 

END M 94- 
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ESTIMATE CHANGES IN INITIAL VALUE GUESSES * N 2 

N A 

DELH(J) = HPARt J , K ) *HF IN ( K , L ) *UPDAT ( L ) N 5 

N 6 

WHERE HPAR = A MATRIX OF PARTIAL DERIVATIVES OF THE TERMINAL N 7 

CONSTRAINTS H(J), ViRT. THE STATES CON STR . ( N-M+ 1 X KS) N 8 

KS = NUMBER OF VARIABLES INVOLVED IN THE CONSTRAINTS. N 9 

N 10 

HFIN = A MATRIX WHOSE FIRST N-M COLUMNS CONTAIN THE KS N 11 

FINAL VALUES OF THE CONSTRAINED VARIABLES RESULTING N 12 

FROM EACH OF THE HOMOGENEOUS SOLUTIONS. N 13 

THE LAST COLUMN CONTAINS THE DERIVATIVES OF THESE N IA 

VARIABLES EVALUATED AT THE END POINT. (KS X N-M+l) N 15 

UPDAT = A VECTOR OF CHANGES TO BE MADE IN THE INITIAL GUESSES N 16 

DELTA-T FINAL IS THE LAST ELEMENT. (N-M+l) N 17 

DELH = TERMINAL BOUNDARY CONDITION DISSATISFACTIONS. (N-M+l) N 18 

N 19 

I F I NAL= DESCRIBES POSITION OF FINAL CONSTRAINED STATES IN X-VECTOR N 20 

N 21 

SUBROUTINE BOUND ( X , DER , CONST , PAR , UP DAT ) N 22 

REAL X(IOO) ,DER( 100) .CONST! 10) ,HPAR( 10, 10) ,HFIN! 10,10) ,H( 10, 10) ,PA N 23 

1R ( 20 ) N 2 A 

REAL DELHI 10 ), UPDAT ( 10 ) N 25 

INTEGER I F INAL ( 10 ! N 26 

REAL LU,LG,LR,M N 27 

N=PAR(2) N 28 

M= P AR ( 3 ) N 29 

I TF=PAR I 5 ) N 30 

NTC-N-M+ ITF N 31 

DO 1 1=1,10 N 32 

DO 1 J=l,iO N 33 

HPARt I , J 1=0.0 N 3A 

N 35 

THE USER MUST PROVIDE SECTIONS TO EVALUATE DE LH, HPAR , I F I NAL . N 36 

THESE MUST DE EVALUATED USING FINAL CONDITIONS. N 37 

FORM TERMINAL CONSTRAINT DISSATISFACTIONS. N 38 

VARIABLES INVOLVED IN CONSTRAINTS HAVE A 1 IN THE CORRESPOND- N 39 

ING POS IT I ON OF IF INAL N AO 

M=X (25 ) N A1 

U=X(l) N A2 

G= X ( 2 ) N A3 

R=X ( 3 ) N AA 

LU=X ( A ) N A 5 

LG= X ( 5 ) N A6 

LR= X ( 6 ) N A 7 

DELH(l)=CCNST(l)-!U**2/2.0-1.0/R) N A8 

DELH(2)=C0NST(2)-(R*U*C0SIG> ) N A9 

DELH(3)=C0NST(3)-( (U/R-1.0/(U*R**2) ) *LG*COS ( G ) + ( LR «U-LU/R**2 ) *S I NT N 50 

1G ) ) N 51 

DELH(A)=CCNST(A)-(LU**2+(LG/U)**2-M**2) N 52 

N 53 

FORM HPARt J, I ) N 5A 

HPAR ( 1 , 1 ) =U N 55 

HPARtl ,3)=1.0/R**2 N 56 

HPAR!2,1)=R*C05(G) N 57 

HPAR(2,2)=-R*U*SIN(G) N 58 

HPAR(2,3)=U*C0S(G) N 59 

HPAR(3, 1 )=( l.O/R+l.O/ ( (U*R)**2) )*LG*COS' G)+LR*SIN(G) N 60 

HPAR(3,2)=-(U/R-1.0/ (U*R**2) ! *L G*S I N ( G ) + ( LR*U-LU/R**2 ) *COS ( G ) N 61 

HPAR(3,3)=(-U/R**2+2./ (U*R**3! )*LG*C0S(G)+2.*LUTSIN(G)/R*A3 N 62 

HPAR(3,A)=— S IN(G)/R**2 N 63 

HPAR(3,5)=(U/R~1.0/(U*R**2))*C0S(G) N 6A 

HPAR(3,6)=U*SIN(G) N 65 

HPARIA, 1 )=-2.0*LG**2/U**3 N 66 

HPARt A, A)=2.0*LU N 67 

HPAR(A,5)=2.0*LG/U**2 N 68 

N 69 

FORM THE DESCRIPTION OF POSITION OF CONSTRAINT VARIABLES. N 70 
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IF INAL 1 1 ) = 1 N 71 

I F I NAL ( 2 ) = 1 N 72 

I F I NAL ( 3 ) = 1 N 73 

IF INAL l A ) = 1 N 76 

IF I NAL 1 5 ) = 1 N 75 

1 FI NAL ( 6 J — 1 N 76 

C N 77 

C EVALUATE THE FINAL MATRIX. N 78 

LL=N-M N 79 

DO 2 I = 1 . LL N 80 

L=0 N 81 

DO 2 K=1,N N 82 

IF ( I F I NAL ( K ) . E Q . 0 ) GO TO 2 N 83 

L=L+1 N 86 

HFINIL. I )=X( I*N+K) N 85 

2 CONTINUE N 86 

C N 87 

IF (ITF.EQ.O) GO TO 4 N 88 

L=0 N 89 

DO 3 K= 1 » N N 90 

IF (IFINAL(K).EQ.O) GO TO 3 N 91 

L=L + 1 N 92 

HFINIL, NTC)=DER(K) N 93 

CONTINUE N 96 

N 95 

FORM PRODUCT N 96 

DO 6 J= 1 , NTC N 97 

DO 6 K= 1 * NT C N 98 

SUM=0. 0 N 99 

DO 5 1=1, L N 100 

SUM= SUM + HPAR (J,I )*HFIN(I,K) N 101 

H ( J , K ) = SUM N 102 

N 103 

SOLVE LINEAR EQUATIONS FDR UPDATES OF INITIAL CONDITIONS. N 106 

UPDAT = HI I NVERSE ) *DELH. (HPAR = DUMMY TO TAKE INVERSE.) N 105 

CALL SIMEQ ( H, DELH, NTC , HP AR .UPDAT ) N 106 

RETURN N 107 

END N 108- 

ESTIMATE CHANGES IN INITIAL VALUE GUESSES * 02 

Q 3 

0 A 

DELH(J) = HPARl J , K ) *HF IN ( K, L ) *UPD AT l L ) 0 5 

0 fc 

WHERE HPAR = A MATRIX OF PARTIAL DERIVATIVES OF THE TERMINAL 0 7 

CONSTRAINTS H(J), WRT. THE STATES CONSTR . IN-M+ 1 X KS) 0 8 

KS = NUMBER OF VARIABLES INVOLVED IN THE CONSTRAINTS. 0 9 

' 0 10 

HFIN = A MATRIX WHOSE FIRST N-M COLUMNS CONTAIN THE KS 0 11 

FINAL VALUES OF THE CONSTRAINED VARIABLES RESULTING 0 12 

FRCM EACH OF THE HOMOGENEOUS SOLUTIONS. 0 13 

THE LAST COLUMN CONTAINS THE DERIVATIVES OF THESE 0 16 

VARIABLES EVALUATED AT THE END PQINT. (KS X N-M+l) 0 15 

UPDAT= A VECTOR OF CHANGES TO BE MADE IN THE INITIAL GUESSES 0 16 

DELTA-T FINAL IS THE LAST ELEMENT. (N-M+l) 0 17 

DELH = TERMINAL BOUNDARY CONDITION DISSATISFACTIONS. (N-M+l) 0 18 

0 19 

IFINAL= DESCRIBES POSITION OF FINAL CONSTRAINED STATES IN X-VECTOR 0 20 

0 21 

SUBROUTINE FIXED ( X , DER, CONST , PAR , UPDAT ) 0 22 

REAL XI 100) ,DERI 100) .CONST! 10) ,HPAR( 10, 10) ,HFIN( 10, 10) ,H(10,10) ,PA 0 23 

1R120) 0 2A 

REAL DELHI10 ) , UPDAT ( 10) 0 25 

INTEGER IF I NAL ( LO J 0 26 

REAL LR , LU , LG , M 0 27 

N=P AR t 2 ) 0 28 

M=P AR ( 3 ) 0 29 

I TF = PAR ( 5 ) 0 30 

NTC=N-M+1TF 0 31 

DO 1 1=1,10 0 32 

DO 1 J=l,10 0 33 

1 HPAR ( I , J ) =0.0 0 36 

C 0 35 
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THE USER MUST PROVIDE SECTIONS TO EVALUATE DELHiHPAR, [FINAL. 0 36 

THESE MUST BE EVALUATED USING FINAL CONDITIONS. 0 37 

FORM TERMINAL CONSTRAINT DISSATISFACTIONS. 0 38 

VARIABLES INVOLVED IN CONSTRAINTS HAVE A L IN THE CORRESPOND- O' 39 
ING POSITION OF I F I HAL 0 AO 

u=xm o 4i 

G=X (2 ! 0 42 

R=X(3) 0 A3 

LU=X(4i 0 AA 

LG=X ( 5 ) 0 A 5 

LR= X ( 6 ) 0 A 6 

M=X (25 ) 0 A 7 

F = P AR ( 20 ) 0 A8 

SQRTD=SQRT( ( LG/U ) **2+LU**2 ) 0 A9 

DELHtl !=CCNST I 1 )-U 0 50 

DELHI2 )=C0NST (2 )-G 0 51 

DELH(3)=C0NST(3)-R 0 52 

DELH14)=C0NSTI4)-(F*SCRTD/M+CU/R-1.0/(U*R**2 ) ) *LG*C0S ( G ) + I LR*U- L U/ 0 53 

1R**2)*SIN(G)-F) 0 5A 

0 55 

FORM HPARI J, I ) 0 56 

HPARt 1 , 1 )=1.0 0 57 

HPAR ( 2 , 2 ) = 1 . 0 0 58 

HP AR ( 3 , 3 ) = l . 0 0 59 

HPAR (A , 1 )=-F*LG**2/ t SORT D*M*U**3 ) + 1 1 .0/R-H.O/t (U*R !**2) ) *LG*COSIG) 0 60 

1 +LR*S I N ( G ) 0 61 

HPAR{4,2)=-fU/R-1.0/IU*R**2> K-LG* S I N ( G ) + ( LR *U-LU/R **2 ) #COS ( G ) 0 62 

HPAR(4,3)=-{U/R**2-2.0/(U*R**3))*LG$C0S(G)+2.0*LU*SIN(G)/R**3 0 63 

HPAR{4,4)=F*LU/(M*SQRTD)-SINtG}/R**2 0 6A 

HP AR 1 4 , 5 ) = F * LG/ ( M*S QRT 0*0 * *2 ) -COS ( G ) * { 1 • 0/ ( U*R** 2 ) ~U/ R ) 0 65 

HPAR(4,6)=U*SIN{G) 0 66 

0 67 

FORM THE DESCRIPTION Of POSITION GF CONSTRAINT VARIABLES. U 68 

IF 1 NAL t 1 1 = 1 0 69 

I F I NAL ( 2 ) = L 0 70 

I F I N AL t 3 ) = 1 0 71 

IF I NAL t A 1 = I 0 72 

I F 1 NAL t 5 ) = 1 0 73 

I F I NAL (63=1 0 7 A 

0 75 

EVALUATE THE FINAL MATRIX. 0 76 

Ll=N-M 0 77 

DO 2 1 = 1, LL 0 78 

L = 0 0 79 

DO 2 K=1,N 0 80 

IF (1FINALIK1.EQ.0) GG TO 2 0 81 

L=L+ 1 0 82 

HFINCL, I )=X( I*N+K) 0 83 

CONTINUE 0 eA 

0 85 

IF (ITF.EQ.O) GO TO A 0 86 

L=0 0 87 

DO 3 K=1 ,N 0 88 

IF (IFINALtK J.EQ.O) GO TO 3 0 89 

L=L+ 1 0 90 

HFIN(L,NTC)=DER(K) 0 91 

CONTINUE 0 92 

0 93 

FORM PRODUCT 0 9A 

DO 6 J= 1 , NTC 0 95 

DO 6 K= 1 i NT C 0 96 

SUM=0. 0 0 97 

DO 5 1=1, L 0 98 

SUM=SUH+HPAR(J, I ) *HF IN ( I , K ) 0 99 

H ( J , K ) = SUM . 0 100 

0 101 

SOLVE LINEAR ECUAT1CNS FOR UPDATES OF INITIAL CONDITIONS. 0 102 

UPDAT = H( INVERSEI-DELH. I HPAR = DUMMY TO TAKE INVERSE.) 0 103 

CALL SIMEQ IHjDELH, NTC, HPAR, UPDAT) 0 10A 

RETURN 0 105 

END 0 106- 
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SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS * P 2 

#**##❖#**#**#**#££#**#❖$##*****:*$$£**&##* P 3 

P A 

XDOT(I) = A(I,J)*X( J) P 5 

P 6 

XDOT(I) = A VECTOR OF KC CONSTANTS. P 7 

KC = THE NUMBER OF LINEAR EQUATIONS. P 8 

A ( I » J ) = A KC BY KC MATRIX OF CONSTANTS. P 9 

SIMEQ SOLVES FOR AND RETURNS THROUGH THE CALLING LIST, P 10 

A I N V = A KC BY KC MATRIX = A! INVERSE) P II 

X(J) = A SOLUTION VECTOR. P 12 

P 13 
P 1 A 

SUBROUTINE SIMEQ ( A , Y , KC , A I NV , X ) P 15 

DIMENSION AIIO.IO), BI10.10), AINV(10,L0), YI10), X ( l 0 ) P 16 

P 17 

SET INVERSE- TO IDENTITY. SAVE A AND Y (THESE ARNT DESTROYED) P 18 

DO 1 1 = 1, KC P 19 

X ( I ) = Y ( I } P 20 

DO 1 J= 1 t KC P 21 

B ( I , J) = A( I , J ) • P 22 

AINV ( I , J )=0. 0 P 23 

IF (I.EQ.J) AINVl I, J)=1.0 P 2A 

CONTINUE P 25 

P 26 

GENERATE INVERSE AND SOLUTION SIMULTANEOUSLY BY TRANSFORMING P 27 

A INTO IDENTITY AND PERFORMING THE SAME OPERATIONS ON I DENT . P 28 

DO 11 1=1, KC P 29 

P 30 

FIND THE LARGEST ELEMENT IN I-TH COLUMN P 31 

C0MP=0. 0 P 32 

DO 3 K= I , KC P 33 

TEMP=B ( K, I ) P 3A 

IF (ABS(TEMP).GT.ABS(COMP) ) 2,3 P 35 

COMP=TEMP P 36 

N=K P 37 

CONTINUE P 38 

P 39 

IF LARGEST ELEMENT IS ZERO, THEN MATRIX IS SINGULAR. P AO 

IF (COMP. EQ. 0.0) A, 6 P A1 

WRITE (3,5) KC , KC P A2 

5 FORMAT!//// 20X, ^LINEAR SIMULTANEOUS EQUATIONS SUESR. SINGULAR P A3 

2 *, 12, * BY *, 12, * MATRIX.*//) P AA 

KC=-KC P A5 

RETURN P A6 

P A 7 

CHECK FOR THE LARGEST ELEMENT ON THE DIAGONAL P A8 

IF NOT ON DIAGONAL, INTERCHANGE COLUMNS I AND N. P A9 

IF (N.EQ. I ) GO TO 8 P 50 

TEMP= XII) P 51 

XII ) = X I N ) P 52 

X t N ) =T E MP P 53 

DO 7 M=1,KC P 5A 

TEMP=8(I,M) P 55 

B(I,M)=8(N,M) P 56 

BIN, M) =TEMP P 57 

TEMP=A I NV ( I , M) P 58 

A I N V 1 1 , M ) = A I N V ( N , M ) P 59 

AINV(N,M)=TEMP P 60 

P 61 
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LARGEST ELEMENT IS NOW ON THE DIAGONAL. P 62 

DIVID THROUGH COLUMN BY THE DIAGONAL ELEMENT. P 63 

TEMP=1 .0/811 » I ) P 64 

xm=X(I)*TEHP P 65 

DO 9 M=1,KC P 66 

Btl ,H:=BtI »P>*TEMP p 67 

AINVI I ,M)=AINV( I rM) *TEMP P 68 

P 69 

DIAGONALIZE B THUS GENERATING AINV AND X. P 70 

DO 11 J=1,KC P 71 

TEMP=BtJ,I) P 72 

IF U.EQ. J.OR.TEPP.EQ.O.O) GO TO 11 P 73 

X(J)=X(J)-TEMP*X(I) P 74 

DO 10 N=l,KC P 75 

8(J*N)=B(J»N)-TEMP*8(I,N) P 76 

10 AINVI J,N)=AINV( J,N)-TEMP*AINV( I,N) P 77 

11 CONTINUE P 78 

RETURN P 79 

END P 80- 
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